1 


AFOSR  FINAL  REPORT 
Grant  No:  FA9550-041-1-0042 

WAYNE  STATE  UNIVERSITY 
Department  of  Mechanical  Engineering 
Detroit,  MI  48202 

Pis:  Raouf  A.  Ibrahim*  and  Ronald  F.  Gibson 
*  Tel.  (313)  577-3885 


Graduate  Students  Supported: 

1.  Dumitru  M.  Beliou 

2.  Stefan  C.  Castravete 

3.  Srinivasa  D.  Thoppul 


Project  Title: 

Nonlinear  Stochastic  Flutter  of  a  Cantilever  Wing  with  Joint  Relaxation 

and  Random  Loading 


Date:  February  21,  2008 


20080331083 


0 


REPORT  DOCUMENTATION  PAGE 

Form  Approved 

OMB  No  0704-0188 

Public  reporting  burden  for  this  collection  of  Information  is  estimated  to  average  1  hour  per  response  including  the  time  for  reviewing  instructions,  searching  eaistng  data 
sources,  gathering  and  maintaining  the  data  needed,  and  completing  and  reviewing  this  oollectlon  of  Information  Send  comments  regarding  this  burden  estimate  or  any 
other  aspect  of  this  collection  of  Information.  Including  suggestions  for  reducing  this  burden  to  Department  of  Defense.  Washington  Headquarters  Services  Directorsts  for 
Information  Operations  and  Reports  (0704-0180).  1215  Jefferson  Davis  Highway,  Suita  1204.  Arlington.  VA  22202-4302  Respondents  should  be  aware  that 
notwithstanding  any  other  provision  of  law,  no  parson  shall  be  subject  to  any  penalty  for  falling  to  comply  with  a  collection  of  information  if  it  does  not  displsy  a  currently  valid 
OMB  control  number  PLEASE  DO  NOT  RETURN  YOUR  FORM  TO  THE  ABOVE  ADORESS. 

1 .  REPORT  DATE  (09-02-2008) 

2.  REPORT  TYPE:  Final 

3.  DATES  COVERED  (From  02-01-04 
-  To  08-30-07) 

4.  TITLE  AND  SUBTITLE 

Nonlinear  Stochastic  Flutter  of  a  Cantilever  Wing  with  Joint  Relaxation  and 
Random  I.oading 

6a.  CONTRACT  NUMBER 

6b.  GRANT  NUMBER 

FA9660-04-1  -0042 

6o.  PROGRAM  ELEMENT  NUMBER 

8.  AUTHOR(S) 

Raouf  A.  Ibrahim,  and  Ronald  F.  Gibson 

Ronald  G.  Gibson 

6d.  PROJECT  NUMBER 

6a.  TASK  NUMBER 

6f.  WORK  UNIT  NUMBER 

7.  PERFORMING  ORGANIZATION  NAME(S)  AND  ADDRESS(ES) 

AND  ADDRESS(ES) 

Wayne  State  University, 

Department  of  Mechanical 

Hngineering,  Detroit,  MI  48202 

8.  PERFORMING  ORGANIZATION 
REPORT 

NUMBER 

9.  SPONSORING  /  MONITORING  AGENCY  NAME(8)  AND  ADDRE8S(ES) 


US.  Air  Force  (AFOSR) 

O  ‘£5>itr 

A/tuvy  •  ~ 


10.  SPONSOR/MONITOR  S 
ACRONYM  (8) 


11.  SPONSOR/MONITOR’ 8 
REPORT 


KA~~  afrl-sr-ar-tr-°8-°151 


12.  DISTRIBUTION  /  AVAILABILITY  STATEMENT: 

Unlimited 


13.  SUPPLEMENTARY  NOTES 

See  the  attached  technical  Report 


14  AB8TRACT 

The  research  project  has  generated  a  number  of  problems  addressing  the  uncertainties  and 
relaxation  problems  in  aeroelastic  structures.  Specifically,  three  main  problems  that  are  of 
important  concern  to  the  aerospace  industry  and  the  Air  Force  technology  have  been 
addressed.  These  are: 

1 .  The  influence  of  structure  uncertainties  of  the  flutter  of  an  aircrafl  wing. 

2.  Influence  of  joint  relaxation  on  the  flutter  of  aeroelastic  structures  (experimental  and 
analytical  investigations). 

3.  Stabilization  of  wing  flutter  via  parametric  excitation. 

1'hese  problems  are  addressed  in  more  details  in  the  attached  technical  report. 


16.  SUBJECT  TERMS 


16.  SECURITY  CLASSIFICATION  OF:  Unclassified 


a.  REPORT 


b.  ABSTRACT 


o.  THIS 
PAGE 


17. 

LIMITATION 
OF  ABSTRACT 


18. 

NUMBER 

OF 

PAGES 


19a.  NAME  OF  RESPONSIBLE 
PERSON 


19b.  TELEPHONE  NUMBER 

(in&udm  mrmm  cod*) 


Standard  Form  298  (Rav.  8- 

98) 

Prescribed  by  ANSI  Sid  Z38.1S 


2 


Table  of  Contents 

REPORT  DOCUMENTATION  PAGE 
AFOSR  FINAL  REPORT 

ABSTRACT . 1 

CHAPTER  1:  EFFECT  OF  STIFFNESS  UNCERTAINTIES  ON  THE  FLUTTER  OF  A 
CANTILEVER  WING . 2 

CHAPTER  2:  THE  INFLUENCE  OF  JOINT  RELAXATION  ON  THE  FLUTTER  OF 

AEROELASTIC  STRUCTURES . 25 

Part  1:  Experimental  Investigation . 25 

Part  2:  Analytical  Investigation . 53 

CHAPTER  3:  FLUTTER  SUPPRESSION  OF  A  PLATE-LIKE  WING  VIA 

PARAMETRIC  EXCITATION . 94 

Appendix:  List  of  Publications  Generated  from  the  Grant .  148 

List  of  Ph.D.  Theses  Generated .  149 


1 


ABSTRACT 

The  research  project  has  generated  a  number  of  problems  addressing  the  uncertainties  and 
relaxation  problems  in  aeroelastic  structures.  Specifically,  three  main  problems  that  are  of 
important  concern  to  the  aerospace  industry  and  the  Air  Force  technology  have  been  addressed. 
These  are: 

1 .  The  influence  of  structure  uncertainties  of  the  flutter  of  an  aircraft  wing. 

2.  Stabilization  of  wing  flutter  via  parametric  excitation. 

3.  Influence  of  joint  relaxation  on  the  flutter  of  aeroelastic  structures. 

These  problems  are  addressed  in  the  next  three  chapters. 
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CHAPTER  1. 

EFFECT  OF  STIFFNESS  UNCERTAINTIES  ON  THE 
FLUTTER  OF  A  CANTILEVER  WING 


Summary  of  Results 

This  research  project  deals  with  the  influence  of  structural  uncertainties  and  joint  relaxation  on  the  flutter 
probability  and  safety  of  nonlinear  cantilever  wings.  The  tasks  of  the  research  have  been  accomplished 
using  analytical,  numerical,  and  experimental  techniques.  The  influence  of  span-wise  distribution  of 
bending  and  torsion  stiffness  uncertainties  on  the  flutter  behavior  of  an  aeroelastic  wing  using  a 
stochastic  finite  element  approach  was  studied  using  a  numerical  algorithm  to  simulate  unsteady, 
nonlinear,  incompressible  flow  (based  on  the  unsteady  vortex  lattice  method)  interacting  with 
linear  aeroelastic  structure  in  the  absence  of  uncertainties.  The  air  flow  and  wing  structure  were 
treated  as  elements  of  a  single  dynamical  system.  Parameter  uncertainties  were  represented  by  a 
truncated  Karhunen-Love  expansion.  Both  perturbation  technique  and  Monte  Carlo  simulation 
were  used  to  establish  the  boundary  of  stiffness  uncertainty  level  at  which  the  wing  exhibits 
flutter  in  the  form  of  limit-cycle  oscillations  (LCO)  and  above  which  the  wing  experiences 
dynamic  instability.  The  results  of  the  perturbation  approach  were  compared  with  those  predicted 
by  Monte  Carlo  simulation  and  the  comparison  revealed  good  correlation  for  low  values  of 
stiffness  uncertainty  levels.  The  major  findings  are: 

•  As  uncertainty  level  increases,  the  perturbation  method  loses  the  accuracy. 

•  For  the  prediction  of  LCO,  the  perturbation  method  is  very  accurate  for  all  levels  of 
bending  stiffness  uncertainty  examined,  but  the  method  loses  its  accuracy  at  upper  levels 
of  torsion  stiffness  uncertainty. 

•  The  stability  boundary  in  the  flow  speed  versus  stiffness  uncertainty  reveals  the 
appearance  of  LCOs  just  below  the  flutter  speed  boundary.  Further  increase  of 
uncertainty  level  produces  instability. 

•  The  uncertainties  in  torsion  stiffness  induce  a  greater  disturbance  in  the  system.  A  smaller 
level  of  torsion  stiffness  uncertainty  induces  instability  in  the  system. 

State-of-the-Art 

The  presence  of  parameter  uncertainties  in  aeroelastic  structures  adds  a  new  dimension  to  an 
already  complicated  problem.  Parameter  uncertainties  owe  their  origin  to  a  number  of  sources, 
which  include: 

(i)  randomness  in  material  properties  due  to  variations  in  material  composition; 

(ii)  randomness  in  structural  dimensions  due  to  manufacturing  variations  and  thermal  effects; 

(iii) randomness  in  boundary  conditions  due  to  preload  and  relaxation  variations  in 
mechanical  joints;  and, 

(iv) randomness  of  external  excitations. 

Generally,  uncertainty  is  described  as  either  parametric  or  non-parametric.  Parametric 
uncertainty  is  due  to  variability  in  the  value  of  input  parameters,  while  non-parametric 
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uncertainty  includes  all  other  sources,  such  as  modeling  errors,  coarse  finite  element  mesh 
fidelity,  or  un-modeled  nonlinear  effects. 

Parameter  uncertainties  may  cause  sensitivity  and  variability  of  the  response  and  eigenvalues  of 
structural  stochasticity1'3.  The  early  developments  relied  on  Monte  Carlo  simulation  and  later  on 
first-  and  second-order  perturbation  methods  to  compute  second-order  moments  of  structure 
response.  Furthermore,  the  general  sources  of  uncertainty  affecting  the  design  and  testing  of 
aeroelastic  structures  were  discussed.  In  particular,  Pettit3  addressed  a  number  of  applications  of 
uncertainty  quantification  to  various  aeroelastic  problems  such  as  flutter  flight  testing,  prediction 
of  limit-cycle  oscillations  (LCO),  and  design  optimization  with  aeroelastic  constraints.  Different 
computational  methodologies  have  been  employed4-5  to  quantify  the  uncertain  response  of 
aeroelastic  structures  with  parametric  variability.  These  methodologies  include  finite  element 
and  perturbation  methods. 

One  of  the  major  problems  of  incorporating  the  random  field  into  finite  element  analyses  is  to 
deal  with  abstract  spaces  which  have  limited  physical  support6,7.  The  difficulty  involves  the 
treatment  of  random  variables  defined  on  these  abstract  spaces.  Usually  the  problem  is  solved 
through  Monte  Carlo  simulation  or  stochastic  finite  element  methods.  Due  to  the  large  number  of 
samples,  which  require  high  computational  time,  the  Monte  Carlo  simulation  is  used  mainly  to 
verify  other  approaches.  The  perturbation8'11  and  Neumann  expansion10'12  methods  proved 
acceptable  results  for  small  random  variation  in  the  material  properties.  It  was  found  that  these 
methods  are  comparable  in  accuracy,  but  the  most  efficient  solution  procedure  is  the  perturbation 
finite  element  method,  which  requires  a  single  simulation.  However,  perturbation  method 
requires  the  system  uncertainty  to  be  small  enough  to  guarantee  convergence  and  accurate 
results. 

The  perturbation  stochastic  finite  element  method  (SFEM)  has  been  adopted  by  several 
researchers  using  the  Karhunen-Loeve  (K-L)  expansion  to  discretize  the  random  fields  due  to 
structure  mechanical  properties13'15.  Jensen16  considered  an  extension  of  the  deterministic  finite 
element  method  to  the  space  of  random  function.  A  Neumann  dynamic  SFEM  of  vibration  for 
structures  with  stochastic  parameters  under  random  excitation  was  treated  by  Lei  and  Qiu17.  The 
equation  of  motion  was  transformed  into  quasi-static  equilibrium  equation  for  the  solution  of 
displacement  in  time  domain.  The  Neumann  expansion  method  was  applied  to  the  equation  for 
deriving  the  statistical  solution  of  the  dynamic  response  within  the  framework  of  Monte  Carlo 
simulation.  The  K-L  expansion  has  proven  to  be  a  powerful  tool  in  modeling  parameter 
uncertainties  in  the  structural  dynamics  community.  However,  it  has  not  been  utilized  in  the 
randomness  and  variability  of  parameters  in  aeroelastic  structures.  The  present  work  is  an 
attempt  to  employ  the  K-L  expansion  to  discretize  the  span-wise  distribution  of  bending  and 
torsion  stiffness  uncertainties  of  an  aircraft  wing. 

Structural  and  material  uncertainties  have  a  direct  impact  on  the  flutter  characteristics  of 
aeroelastic  structures  and  they  have  begun  to  attract  some  attention  in  the  literature.  They  were 
considered  in  studying  the  flutter  of  panels  and  shells18'22.  Liaw  and  Yang18,19  quantified  the 
effect  of  parameter  uncertainties  on  the  reduction  of  the  structural  reliability  and  stability 
boundaries  of  initially  compressed  laminated  plates  and  shells.  For  buckling  analysis,  the 
uncertainties  were  included  in  the  modulus  of  elasticity,  thickness,  and  fiber  orientation  of 
individual  lamina,  as  well  as  geometric  imperfections.  For  flutter  analysis,  further  uncertainties 
such  as  mass  density,  air  density,  and  in-plane  load  were  also  considered.  Kuttenkeuler  and 
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Ringertz20  performed  an  optimization  study  of  the  onset  of  flutter,  with  respect  to  material  and 
structural  uncertainties  using  finite  element  analysis  and  the  doublet-lattice  method.  Lindsley  et 
al.2  ’  2  considered  uncertainties  in  the  modulus  of  elasticity  and  boundary  conditions  for  a 
nonlinear  panel  in  supersonic  flow.  The  probabilistic  response  distributions  were  obtained  using 
Monte  Carlo  simulation.  It  was  reported  that  uncertainties  have  the  greatest  nonlinear  influence 
on  LCO  amplitude  near  the  deterministic  point  of  LCO.  Poirion23,24  employed  a  first-order 
perturbation  method  to  solve  for  the  probability  of  flutter  given  uncertainty  in  the  structural  mass 
and  stiffness  operators. 

Recently,  the  influence  of  parameter  uncertainties  on  the  response  of  a  typical  airfoil  section  was 
considered  by  few  researchers25'28.  The  sensitivity  of  uncertainty  of  aeroelastic  phenomena  was 
evaluated  using  Monte  Carlo  simulation.  The  effect  of  parametric  uncertainty  on  the  response  of 
a  nonlinear  aeroelastic  system  was  studied  by  Attar  and  Dowell28  using  a  response  surface 
method  to  map  the  random  input  parameters  to  the  root-mean  square  wing  tip  response. 

Civil  engineers  have  been  involved  in  studying  the  influence  of  uncertainties  of  structural 
properties,  in  particular  damping,  on  the  reliability  analysis  of  flutter  of  a  bridge  girder  and  a  flat 
plate  was  determined  in  few  studies29'31.  The  prediction  of  the  flutter  wind  speed  was  found  to  be 
associated  with  a  number  of  uncertainties  such  that  the  critical  wind  speed  can  be  treated  as  a 
stochastic  variable.  The  probability  of  the  bridge  failure  due  to  flutter  was  defined  as  the 
probability  of  the  flutter  speed  exceeding  the  extreme  wind  speed  at  the  bridge  site  for  a  given 
period  of  time.  The  probabilistic  dynamic  response  of  a  wind-excited  structure  has  been  studied 
in  terms  of  uncertain  parameters  such  as  wind  velocity,  lift  and  drag  coefficients  by  Kareem32. 
The  influence  of  uncertainty  in  these  parameters  was  found  to  propagate  in  accordance  with  the 
functional  relationships  that  relate  them  to  the  structural  response.  Note  that  while  in  aerospace 
structures,  one  is  interested  in  estimating  the  onset  flutter  due  to  parameter  uncertainty,  civil 
engineers,  on  the  other  hand,  focus  on  probabilistic  reliability  analysis  to  determine  a  probability 
of  the  bridge  or  a  structure  failure  due  to  flutter  for  a  given  return  period  rather  than  stating  a 
single  critical  wind  speed. 

A  ground  vibration  test  was  used  by  Potter  and  Lind33  to  obtain  uncertainty  models,  such  as 
natural  frequencies  and  their  associated  variations,  which  can  update  analytical  models  for  the 
purpose  of  predicting  robust  flutter  speeds.  Different  norm  approaches  were  used  to  formulate 
uncertainty  models  that  cover  the  entire  range  of  observed  variations.  Lind  and  Brenner34 
introduced  a  tool  referred  to  as  the  “flutterometer”  for  predicting  the  onset  of  flutter  during  a 
flight  test.  The  flutterometer  computes  the  onset  of  flutter  for  an  analytical  model  with  respect  to 
an  uncertainty  description.  Brenner35  considered  a  technique  that  identifies  model  parameters  and 
their  associated  variances  from  flight  data.  Later,  Prazenica,  et  al.36  introduced  a  technique  for 
estimating  uncertainty  descriptions  based  on  a  wavelet  approach,  but  relies  on  Volterra  kernels. 

The  present  work  deals  with  the  influence  of  stiffness  uncertainties  on  the  flutter  behavior  of  an 
aeroelastic  wing.  A  numerical  algorithm  originally  developed  by  Predikman  and  Mook37  to 
simulate  unsteady,  nonlinear,  incompressible  flow  interacting  with  linear  aeroelastic  wing  in  the 
absence  of  uncertainties  is  adopted.  In  order  to  implement  this  algorithm  in  the  presence  of 
uncertainties  we  introduce  a  random  field  that  represents  bending  or  torsion  stiffness  parameters 
or  both  as  a  truncated  Karhunen-Love  (K-L)  expansion6.  The  air  flow  and  wing  structure  are 
treated  as  elements  of  a  single  dynamic  system.  Both  perturbation  technique  and  Monte  Carlo 
simulation  are  used  to  establish  the  boundary  of  stiffness  uncertainty  level  at  which  the  wing 
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exhibits  LCO  and  above  which  the  wing  experiences  dynamic  instability.  In  this  paper  the  term 
limit  cycle  oscillation  is  used  to  denote  the  flutter  boundary,  where  for  a  linear  system,  effective 
damping  is  zero.  The  analysis  also  includes  the  limitation  of  perturbation  solution  for  relatively 
large  level  of  stiffness  uncertainty.  The  analytical  modeling  of  aerodynamic  loading  based  on  the 
unsteady  vortex  lattice  method,  structural  forces  interacting  with  the  aerodynamic  loading,  and 
stiffness  uncertainties  based  on  The  K.-L  expansion  are  briefly  described  in  sections  II  through 
IV,  respectively.  Section  V  establishes  the  entire  system  modeling  by  combining  the  three 
models  (aerodynamic,  structure,  and  uncertainty)  in  the  wing  governing  equations  of  motion  in 
the  finite  element  discretized  form.  Sections  VI  and  VII  present  the  perturbation  analysis  solution 
and  Monte  Carlo  simulation  results,  respectively,  together  with  a  comparison  of  the 
corresponding  results. 


II.  AERODYNAMIC  MODELING 

Aerodynamic  modeling  requires  the  estimation  of  aerodynamic  forces  and  at  the  same  time 
accounts  for  a  wing  elastic  deformation.  This  is  achieved  by  using  the  unsteady  vortex  lattice 
method  (VLM).  The  VLM  accounts  for  aerodynamic  nonlinearities  associated  with  the  angle  of 
attack,  static  deformation,  vorticity-dominated  flow,  and  unsteady  behavior.  It  is  also  not  limited 
to  small  periodic  motion.  The  general  unsteady  vortex-lattice  model  imitates  the  boundary  layers 
and  the  wakes  as  vortex  sheets37.  The  vortex  sheets  are  of  two  types,  bound-vortex  sheet,  and 
free  vortex  sheet.  The  bound-vortex  sheets  create  the  boundary  layer  on  the  surface  of  the  wing 
while  free-vortex  sheets  represent  the  wakes. 

The  lifting  surface  is  approximated  by  a  set  of  lattices  of  short  segments  of  constant  circulation 
each.  Each  segment  occupies  a  portion  of  the  wing  surface  and  is  enclosed  by  a  loop  of  vortex 
segments.  Figure  1(a)  shows  the  discretization  of  the  lifting  surface,  where  N-frame  is  the  ground 
fixed  coordinate  system  and  B-frame  is  the  coordinate  system  on  the  wing  body  and  is  moving 
with  the  wing.  The  leading  segment  of  the  vortex  loop  is  located  at  a  distance  equivalent  to 
quarter  of  the  panel  length.  The  velocities  are  calculated  at  a  finite  number  of  points  called 
control  or  collocation  points  located  at  the  lattice  center.  Figure  1(b)  represents  the  position  of  an 
arbitrary  control  point  P  on  the  lattice  as  a  result  of  the  wing  structure  displacement. 

The  point  P°  represents  the  point  P  before  the  wing  is  deformed.  In  the  N-frame,  the  positions 
of  P°  and  P  are  given  by  the  vectors,  R0,  and  R  ,  respectively.  The  corresponding  vectors  in  the 
B-frame  are  r0  and  r .  The  displacement  of  the  aerodynamic  point  P  due  to  deformation  of  the 
wing  is,  Ar^  .  The  vectors  R  and  r  can  be  represented  as  functions  of  the  wing  displacement, 

Ar^,  i.e.,  R  =  R0  +  Ar^  ,  and  r=  r0  +  Ar^.  The  pressure  p(t)  at  point  P  was  estimated  using 

Bernoulli’s  equation  for  unsteady  flow  from  which  the  non-dimensional  aerodynamic  load  at 
point  P  is38 


F,,(r)=2(4V[v„-«5*  -*F -S' *'  r']  +  ^.[C,]|j,»  (1) 

where  n  is  the  unit  normal  vector  at  the  control  point,  Ap  =  APIL2C  is  the  non-dimensional  area  of 
the  element  P ,  Ap  is  the  area  of  the  element,  Lc  represents  the  chord-wise  length  of  one  element 
on  the  bound  lattice,  x  =  t/Tc,  R  =  R  /  Z,c  =  R0  +  Ar^  ,  R0  =  R0/LC,  Ar^,  =Ar AJLC,  TC  =  LCIV„  is  a 
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characteristic  time,  VF  =  VF/Vxt  VF  is  the  absolute  velocity  of  the  air  flow,  Vx  is  the  free  stream 
velocity.  AV,  =AV,/VX  is  the  non-dimensional  tangential  velocity  difference  across  the  vortex 
lattice  and  \mf  is  the  “mean”  velocity  which  does  not  recognize  the  presence  of  the  local 
vorticity.  VmF  can  be  considered  as  the  flow  velocity  at  the  midpoint  of  the  vortex  sheet 
thickness.  GP  is  the  circulation  loop  in  non-dimensional  form  of  the  vortex  lattice  element  which 
encloses  the  control  point  P ,  wvfl(r)  is  the  absolute  velocity  of  Ob,  Brp  and  *F  are  the  position 
and  velocity  of  P  relative  to  the  B-frame,  respectively,  and  "ei *(r)  is  the  angular  velocity  of  the 
B-frame.  The  expression  for  the  aerodynamic  load  at  point  P  in  dimensional  form  is: 


Figure  1.  (a)  Lifting  Surface  Discretization,  (b)  The  position  of  an  arbitrary  point  P  on  the  lifting  surface  caused  by 

the  wing  structure  deformation 


F,.,(0  =  [(pX2/2)4]f,.,(t) 


(2) 
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where  p  is  the  density  of  the  fluid.  This  aerodynamic  force  is  acting  at  the  center  of  each  lattice 

and  interacts  with  the  elastic  and  inertia  forces  of  the  wing  structure  and  their  modeling  is  given 
in  the  next  section. 


III.  STRUCTURAL  MODELING 

Figure  2(a)  shows  a  cantilever  wing  with  a  straight  elastic  axis  at  distance  S3  from  the  inertia 
axis.  The  wing  is  structurally  modeled  as  an  Euler-Bernoulli  beam.  The  coordinate  system 
coincides  with  the  B-frame  coordinate  system  from  the  aerodynamic  modeling.  The  wing  is 
divided  into  elements  as  shown  in  Figure  2(b),  where  a  single  element  with  two  nodes  is  shown 
in  Figure  2(c). 

The  displacements  of  each  element,  w(.y,/),  \(y,i),  and  \v(y,t),  are  expressed  in  terms  of  the 
nodal  displacements  u, ,  v,,  w,  in  x,  y,  and  r  directions,  respectively,  and  nodal  rotations  /?, , 
a, ,  yl  about  x ,  y ,  and  r  axes,  respectively.  The  relationships  between  the  element  and  nodal 
displacements  are 


«(*0 = Y«  (y)ue  (') .  = Yv  tv)- v* (0 

*M  =  K  (y)we  (0 ,  &(y,t)  =  \Ta  (y)ae  (/) 

where  (y)  =  (>-)  =  (y3  Ys  Yt  r6},  (y)  =  \ra(y)  =  {Y,  y2). 


(3) 


U.  = 


f  ^ 

U, 


Y/ 

Y/+i 


v  = 


w.W- 


p, 


•  -M-fc 


Yj,  7  =  1.2 . 6  the  shape  functions  given  by  Preidikman39 .  Applying  Lagrange’s  equation  to  each 

coordinate  gives  the  following  set  of  equations  of  motion: 


y'*r  y‘*r  /  \f  )  ( y>+1 

*  Jy.ooy.oO>  M0+  j£/fY„oo'(Y.oo')  dy  u,(/)=  \F,\Sy)dy 

n  y<  /  V  y>  )  v  y>  j 


m 


\  y> 

r  y>* i 


iv.MY.W^  v,!^  J^£YvOO'(yvOO')  dy  v.(/)-  "f  FvYy(y)dy 

<  y>  J  \y,  J  \y, 


m 


)  \(y)YAy)Tdy 


\  y> 


w,  (/)  + w83 


y>*  i 

Jy.O0Y.O0> 


v  y> 


)  Wy)Ya(y)r  cfy 


V  yi 


d„(t)  +  md} 


fi+ 1 

Jy.O0Y.O0> 


v  y> 


*.(') 


(4a) 

(4b) 


]  ^Y.OO’fY.OO'  Jdy  w.  (/)  =  f'|  FwYw(y)dy\  (4c) 

\y>  y  \y,  J 


y,+i  /  \r  \  [  ' 

J  cGJYa  (y)'  ( Y0  (y)’  J  dy  «.  (f) =  J  Ma  Y„  (y)  dy  (4d) 

V*  )  \  y,  j 
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Figure  2.  (a)  Schematic  diagram  of  a  cantilever  wing  model,  (b)  Finite  elements  discretization,  (c)  A  finite  element 

representation 

where  a  prime  denotes  a  derivative  with  respect  to  spatial  variable  y  .  m  is  mass  per  unit  length, 
/0  is  mass  moment  of  inertia  (about  inertia  axis)  per  unit  length,  and  l  is  the  wing  length.  A  is 
the  wing  cross-sectional  area,  E  is  Young’s  modulus,  Ix  is  the  area  moment  of  inertia  of  the 
wing  cross-section  about  x  axis,  I,  is  the  area  moment  of  inertia  of  the  wing  cross-section  about 
r  axis,  J  is  the  polar  moment  of  inertia  of  the  wing  cross-section  about  r  -axis,  G  is  the 
modulus  of  rigidity,  and  c  is  a  factor  accounting  for  the  non-circular  geometry  of  the  wing  cross- 
section.  The  value  of  c  depends  on  the  cross-section  geometry  and  is  documented  in  the 
literature40.  Fu ,  Fv,  and  F„  are  the  external  forces  acting  along  x-,  y-,  and  z-axes,  respectively, 
and  Ma  is  the  external  moment  about  the  inertia  axis.  Equations  (4)  may  be  written  in  the  matrix 
form: 

M'U'+K'U' =F'  (5) 

where  Ue  =  {«,  v,  w,  P,  a,  Y/  «j+l  v(+1  w(+l  P(+1  a(tl  Yj+1}r. 

Equation  (5)  constitutes  the  governing  equations  of  motion  for  a  single  element  without 
uncertainties.  The  inclusion  of  material  uncertainties  will  be  considered  in  the  next  section. 
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IV.  MODELING  OF  MATERIAL  UNCERTAINTIES 

Material  uncertainty  is  considered  for  bending  and  torsion  stiffness  parameters.  Let  the  bending 
and  torsion  stiffness  parameters  be  represented  in  the  form, 


EIX  ( y )  =  El x+  El x  ( y ) ,  GJ(y)  =  GJ+GJ(y) 


(6a) 


The  mean  values  E~IX  »  0 ,  and  GJ  »  0  are  assumed  to  be  much  larger  than  the  root-mean-square 
of  the  random  field  variability  represented  by  Elx(y)  and  GJ(y).  Both  EIx(y)  and  GJ(y)  are 
assumed  to  be  Gaussian  distributed  with  zero  mean  and  their  standard  deviations  aEI  ,  and  aGJ 

are  much  smaller  than  the  corresponding  mean  value,  i.e.,  aE/x /Ylx  « 1 ,  and  uGJ  IGJ  « 1 .  This 
implies  that  the  stiffness  parameters  EIx(y)  and  GJ(y )  form  positive-valued  random  fields.  Let 
the  stiffness  parameters  be  represented  by  the  random  function  x(y$),  where  9  is  a  parameter 
that  belongs  to  the  space  of  random  events  and  ye[-L/2,L/2] .  The  random  field  x{y,&)  can  be 
expressed  by  the  truncated  K-L  expansion6: 

x(y,Q)  =  x(yhi^)ylKf,(y)  (6b) 

n= 1 


where  x(y)  is  the  mean  value  of  x(y,9) ,  K  are  the  eigenvalues  of  the  random  field  as  defined 
in  equation  (7)  and  (9b)  below,  /„(y)  is  a  set  of  eigen-functions,  and  £„(G)  is  a  set  of  random 
variables  with  zero  mean  and  E[^  (e)5„  (6)]  =  8^ ,  8„m  is  the  Kronecker  delta.  Where  A„  are  the 
eigenvalues  of  the  covariance  operators  in  equations  (7)  and  (8),  and  /„  (y)  are  the  corresponding 
eigenfunctions  by  the  solution  of  the  integral  equation 

LI  2 

f  C(y,yl)fM)dy,=^fl,(y)  (7) 

-LI  2 


where  C(y,y,)  is  the  covariance  kernel  of  the  random  field  xLv.9)- 

Expansion  (6b)  is  mathematically  well  founded  and  is  guaranteed  to  converge.  The  convergence 
and  accuracy  of  the  K-L  expansion  were  proven  by  Huang  et  al.41  by  comparing  the  second-order 
statistics  of  the  simulated  random  process  with  that  of  the  target  process.  It  was  shown  that  the 
factors  affecting  convergence  are  mainly  the  ratio  of  the  length  of  the  process  over  correlation 
parameter  and  the  form  of  the  covariance  function.  The  K-L  expansion  has  an  advantage  over  the 
spectral  analysis  for  highly  correlated  processes.  For  long  stationary  processes,  the  spectral 
method  is  generally  more  efficient  as  the  K-L  expansion  method  requires  substantial 
computational  effort  to  solve  the  integral  equation.  In  addition,  it  is  optimum  in  the  sense  that  it 
minimizes  the  mean  square  error  resulting  from  truncating  the  series  at  a  finite  number  of  terms. 
The  bending  and  torsion  stiffness  parameters  will  be  modeled  by  one-dimensional  Gaussian 
random  field  models  with  bounded  mean  squares.  Ghanem  and  Spanos6  and  Loeve42  showed  that 
for  a  Gaussian  process  the  K-L  expansion  converges.  The  covariance  kernel  of  the  random  field 
x(y.Q)  may  be  assumed  in  the  form: 
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C(y,y\ )  =  <r2xe^y-y'\"'"  (8) 

where  a]  is  the  variance  of  the  random  field  x  ,  such  that  ax  «  x(y)  implying  that  x(y,&)  will 
always  be  positive,  lmr  is  correlation  length  such  that  lcor  ->  L .  For  one-dimensional  case,  the 
eigenvalue  problem  (7)  possesses  the  closed  form  analytical  solution6 


fn{y)  = 


cos  (a>„y) 

L  sin(2<ynL/2) 
2+  2con 


2 

"  ‘  *1  +1  "lr 


(9a,  b) 


where  con  are  the  roots  of  the  characteristic  equation: 


1  t  (  L) 

- 0)  tan  co— 

1  (  L\ 

(o  H - tan  co — 

[lcor  V  2 )_ 

[  hor  l  2JJ 

(10) 


Introducing  the  expressions  (9)  into  (6)  the  random  field  takes  the  form: 


xUe)  =  xW+S2^(e)°x 


l  l cor 


«= 1 


|K  +  l//L)[^„+sin(2conZ,/2)] 


cos(wn^) 


(11) 


For  bending  stiffness,  the  random  field  x  is  denoted  by  EIx(y )  with  mean  value  EIX  and 
variance  <r\,x ;  and  for  torsion  stiffness  the  random  field  x  is  denoted  by  GJ(y)  with  mean  value 

GJ  and  variance  a2GJ .  For  simplicity,  it  is  assumed  that  both  parameters  are  uncorrelated.  The 
elemental  stiffness  expression  Ke  is  , 


K'  (0)  =  K‘0  +  £  K^„  (0)+  £  K^„  (0) 


where 


*^*♦1  /  ^‘*2  /  \T  _ <?«•!  /  \T  _ X+i 

ks  =  £/,  }y;(y;)  Jy;(y;)  <*<+£/,  |y/(y„*)  dy  +cgj  }  y;y0'7^ 

y, 


1 

/ 


(12) 

(13a) 


Jki 


*M  =  j2<f„(*K, 


(oJL. 


\/lir)[Lco„+sin(2(onL/2)] 


COS 


(MY. 


’Y'7c£ 

t  w  s 


(13b) 


•V-i 


K‘„=  \2^n(0)aGJ 


coJL. 


(co2n  +l/l2or)[Lo)„  +sin(2<y„£/2)] 


cos(^)Y0'Ya'7^ 


(13c) 


Assembling  the  elemental  matrices  we  obtain 
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MU(0)  + 


■ 

(  N  > 

(  N  y 

K0  + 

2X.U8) 

4- 

IX«U0> 

H 

1  J 

\n- 1  J  _ 

u(e)  =  F,(u) 


(14) 


where  M  is  a  (6 ns  x  6 ns )  mass  matrix  represents  the  assembled  elemental  mass  matrix  M" ,  ns 
is  the  number  of  points  in  the  structural  grid,  K0  is  a  (6«sx6«s)  matrix  representing  the 
assembled  elemental  mean  value  of  the  stiffness  matrix  Kg.  KA(I  and  K,„  are  (6nsx6ns) 
matrices  representing  the  assembled  elemental  bending  Keb„  and  torsion  K‘n  random  stiffness 
matrices,  respectively.  F,  (U,/,0)  is  a  (6wsxi)  vector  representing  the  assembled  elemental  force 
vector  F/ ,  and  U(/,0)  is  a  (6ns  xi)  vector  of  the  displacements  of  the  points  in  the  structural  grid 
and  represents  the  assembled  elemental  displacement  vector  U* . 

Note  that  the  left-hand  side  of  equation  (14)  constitutes  a  set  of  stochastic  second-order 
differential  equations  with  random  variable  coefficients  with  Gaussian  distribution.  In  the 
absence  of  aerodynamic  forces,  these  equations  are  always  stable,  since  the  stiffness  matrices  are 
real  positive  definite.  This  is  guaranteed  by  the  fact  that  the  uncertain  components  of  the  stiffness 
parameters  are  very  small  compared  with  their  mean  values  as  stated  in  the  beginning  of  this 
section.  In  the  next  section  the  aerodynamic  nodes  will  be  connected  to  the  structural  finite 
element  mesh  and  the  aerodynamic  forces  will  be  transferred  to  the  structural  grid. 

V.  SYSTEM  MODELING 

Figure  3(a)  shows  the  structural  grid  superimposed  to  the  aerodynamic  grid. 

The  relationship  between  the  displacements  of  the  aerodynamic  grid  points  and  displacements  of 
the  structural  grid  points  is: 


**A 


Ar, 

Ar2 


=  G^U 


Ar. 


nA  J 


(15) 


where  A9?a  is  a  (3/j^xl)  vector  representing  the  displacements  of  the  control  points  in  the 
aerodynamic  grid,  nA  is  the  number  of  control  points  in  the  aerodynamic  grid,  GAS  is  an  3nA  x  6ns 
interpolation  matrix  that  connects  the  displacements  of  the  nodal  points  in  the  aerodynamic  grid 
to  the  displacements  of  the  nodal  points  in  the  structural  grid37. 

Figure  3(b)  shows  how  an  aerodynamic  point  is  connected  to  an  internal  structural  grid  point, 
where  an  aerodynamic  point  PA  is  connected  to  a  point  Ps  on  the  elastic  axis  between  two 
successive  structural  nodes.  For  the  present  study,  the  structural  nodes  are  labeled  by  the  index  / , 
where  /  =  l,2,...,8(=ws)  while  the  aerodynamic  points  are  labeled  by  the  index  j,  where 
7  =  1,2,...,18(=«a).  The  points  Ps  and  PA  are  in  the  same  plane,  which  is  perpendicular  to  the 
elastic  axis.  The  relative  position  between  PA  and  Ps  is  given  by  the  vector  rfs  =[xj,yj,zJY  The 

displacement  of  the  point  PA  is  given  by  Ar/  ^Ar/.Ar/.A//}7  and  the  displacement  of  the  point 
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Ps  is  given  by  U,  =  a*,  .  The  relationship  between  the  points  PA  and  Ps  is  given 

by38: 


Structural  Nodes 


Aerodynamic  Nodes 

Structure  nodes:  i  =  l,2,...,8,  Aerodynamic  control  points:  j  =  l,2,...,18 

(a) 


Figure  3.  (a)  Structural  grid  superimposed  to  the  aerodynamic  grid,  (b)  Connection  of  an  aerodynamic  grid  point  to 

an  internal  point  of  the  structural  grid  point 


=C„AJ, 


1  0  0 
0  1  0 
0  0  1 


0 

■z 

0 


0 


-x, 


0 

XJ 

0 


(16) 


The  displacement  vector  is  obtained  by  finite  element  interpolation  as  a  function  of 
displacements  on  nodes  i  and  /  + 1  as  follows37: 


U,  = 


■12 


*13 


.Yll  ^22  ^23  ^24. 


«W=YUI/+1 


(17) 


13 


Where  UM+1  =  {«,  v,  w,  p,  a,  y,  w,+1  vj+1  wi+l  p,+l  a,+l  y(+1}7  ,  and  the  detailed  structures 
of  Y,y  are  given  in  the  Appendix.  Introducing  equation  (17)  into  equation  (16)  gives: 


Next,  the  matrix  GJt  will  be  assembled  into  the  global  matrix,  G^. ,  which  is  introduced  into 
equation  (15).  In  order  to  obtain  the  relationship  between  the  aerodynamic  forces  F<  and 
structural  forces  Fs,  the  two  force  systems  must  have  the  same  work  for  any  virtual 
displacement,  i.e.. 


5^  )T  Fa  =  5U£  {(gas f  FA }  =  8Uj,Fs  ( 1 9) 

where  8(A9t„)  and  8u£,  are  virtual  displacements  of  control  points  of  the  aerodynamic  grid  and 
of  structural  nodes,  respectively.  G^  is  the  global  matrix  which  connects  the  control  points  of 
the  aerodynamic  grid  to  the  structural  grid  points,  F^  is  the  matrix  of  aerodynamic  force,  which 
is  located  at  the  aerodynamic  control  point,  and  Fs  is  the  matrix  of  the  structural  load,  which  is 

located  at  the  structural  grid  points.  Equation  (19)  yields  the  relationship  between  structural 
forces  and  the  aerodynamic  forces, 

Fs=(Gy4SfF^=G^  (20) 

where  (g^  =  Gsa  .  Introducing  equation  (20)  into  the  equation  (14),  the  equations  of  motion  of 
the  system  is: 


Mi)(e)+ 

■ 

(  N 

(  N 

y 

K0  + 

ZK/,.„Ue) 

+ 

2X„Ue) 

_ 

1 

) 

Kn= 1 

)  _ 

u(e)=c^F,(e) 


(21) 


Equation  (21)  constitutes  the  equations  of  motion  of  a  wing  involving  stiffness  uncertainties.  A 
modal  analysis  of  equation  (21)  is  performed  and  gives: 

K0V  =  MVA  (22) 

where  T  is  an  N  x  N  eigenvector  matrix  where  N  is  the  number  of  degrees  of  freedom  of  the 
finite  element  model,  and  A  is  a  NxN  diagonal  matrix  of  eigenvalues.  Equation  (21)  will  be 
solved  using  perturbation  analysis  in  the  next  section. 

VI.  PERTURBATION  ANALYSIS 

In  order  to  develop  the  perturbation  analysis,  it  is  convenient  to  introduce  the  following 
coordinate  transformation,  which  approximates  the  deformation  of  the  beam: 
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N 


UM)  =  2>*<7*='Pq; 

k= 1 


(23a) 


where  T  is  a  truncated  form  of  Y  .  For  the  present  application,  the  method  requires  displacement 
vector  be  expanded  into  the  truncated  Taylor  series: 


U(f,0)  =  »P 


*io(')+Z|rlt=o  Z^7Trk=o  «,(*>&,(*)+•• 

n=\  °$n 


2 a 


(23b) 


Note  that  V  is  an  Nxiv  matrix  that  includes  first  N  mode  shape  vector,  q  is  an  N* l 
generalized  coordinates  vector,  and  q0  is  the  mean  value  of  q.  Substituting  equation  (23b)  into 
equation  (21),  pre-multiplying  both  sides  by  M"14'T  and  collecting  terms  of  the  same  power  of  S, 
yields  a  set  of  perturbational  equations  of  zero-order,  first-order,  second-order,  and  so  on.  For 
simplification  sake  we  restrict  the  analysis  up  to  the  first-order  and  the  following  perturbational 
equations  are  obtained: 

•  Zero-order  equation  (coefficients  of  ): 

q0+[vrK0‘F]q0=[M-,'FTG^]F/(-0  (24a) 

•  First-order  equation  (coefficients  of  ): 


q,n  +[*K7  K0*F  Jq,_„  =  -[*TKM*  +  'FTK,,„'K]q0  +  [m-,'I'TG&4]|^. 

J  L  J  04 „ 


(24b) 


where  n  =  \...N,  qu=^L,  and  j>rK0'l']  is  am  NxN  diagonal  matrix  of  first  N  eigenvectors. 

The  aerodynamic  loading  as  defined  by  equation  (1)  depends  on  the  wing  displacement  Ar^ 

which  is  time  dependent.  Note  that  ^4.  js  obtained  by  taking  the  derivative  of  equation  (1)  with 

respect  to  after  substituting  for  Ar^  which  is  expressed  in  equation  (16)  and  equation  (23b) 
which  is  a  function  of  . 


Equations  (24)  are  numerically  integrated  once  using  an  adapted  Hamming’s  fourth  order 
predictor-corrector  method.  The  numerical  solutions  are  performed  using  the  following  values  of 
the  wing43:  EFx=l06Nm2;  EIZ  =50xl06Mw2 ;  g7  =  1.5x106M«2;  EA  =  20x\06N;  m=  10  kg/m' 
l0=\5kgm;  83  =  0.15/w;  wing  length  L  =  3m;  wing  chord  C  =  lm;  angle  of  attack  afl=5°; 
Lc  =  0.16667/m;  and  lcor=L  =  3m.  In  order  to  obtain  realistic  results,  the  wing  will  be  discretized 
using  9  elements  for  the  wing  structure  and  6x10  elements  for  the  lattice.  The  response  will  be 
obtained  from  perturbation  method  at  various  flow  speeds  and  various  bending  and  torsion 
stiffness  uncertainty  levels,  i.e., 
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Qben  (r>^)  tfo.ben  (r)  .hen  (r)^n  ’  9W  ~  tfo.lor  (r)  +  ^  4l,lor  (*")^n 

n*l  ns] 


(25a,b) 


Figure  4.  Dependence  of  wing  response  variances  on  the  number  of  terms  N  in  the  Karhunen-Loeve  expansion 

showing  the  convergence  is  achieved  for  N  >2  . 

where  q0be„  and  q0lor  are  obtained  from  the  equation  (24a),  ql  hen  and  qllor  from  the  equation 

(24b).  The  analysis  is  carried  out  for  N  =  2  since  for  N  >  2  the  results  converge  as  shown  in 
Figure  4.  The  temporal  component  of  the  response  variance  along  the  ensemble  is  given  by: 

SL  ('><?)  =  {iben  (i,0)-E[qhen  (t,0)])2 ,  si  {t,6)  =  (q,or  (t,0)~  E[q,or  (r,0)])2  (26a, b) 

The  wing  bending  and  torsion  time  history  records  are  numerically  estimated  using  equations 
(25)  together  with  the  corresponding  variances  (as  computed  from  equations  (26))  for  airflow 
speed  120  m/s  and  angle  of  attack  a  =  5°,  in  the  absence  of  uncertainty.  The  solutions  are 
obtained  for  initial  conditions  ^„(0)  =  -0.05  ,  qhen( 0)  =  0.0,  9(or(0)  =  0.01 ,  and  <7,„r(0)  =  0.0  .  It  is  found 
that  the  wing  is  stable  under  this  airflow  speed.  In  the  presence  of  small  level  of  bending  stiffness 
uncertainty  up  to  <t£/j  /  EIX  =  0.09  and  zero  torsion  uncertainty,  a0J/GJ  =  0.0,  and  under  the  same 

airflow  speed  the  wing  remains  stable  and  both  time  history  records  and  their  variances  decay 
with  time.  When  the  bending  stiffness  uncertainty  reaches  the  critical  value  aKl  IEIX=  0.1  the 

wing  experiences  instability,  loses  its  zero  equilibrium  state  and  reaches  the  flutter  boundary  in 
the  form  of  LCOs  due  to  aerodynamic  nonlinearity.  Above  that  level  of  bending  stiffness 
uncertainty  the  wing  oscillations  grow  without  limit.  As  the  air  flow  speed  increases,  the  flutter 
boundary  takes  place  for  lower  values  of  bending  stiffness  uncertainty  defined  by  the  bound 
I  El x<  0.1 .  An  extensive  number  of  numerical  solutions  have  been  carried  out  for  the  purpose 

of  establishing  the  flutter  boundary  on  air  flow  speed,  Vx ,  versus  the  variance  of  bending 
stiffness  uncertainty,  <r£/  /  EIX.  Figure  5(a)  shows  the  stability  bifurcation  diagram  <te,  /eFx, 
where  the  region  occupied  by  small  empty  circles,  “o”,  designates  stable  wing  response,  and  the 
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region  covered  by  small  empty  triangles,  “A”,  belongs  to  an  unstable  wing.  The  line  separating 
the  two  regions  signals  the  occurrence  of  flutter,  designated  by  empty  squares 
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Figure  5.  (a)  Stability  bifurcation  diagram  Vx  vs.  cth  /  EIX ,  (b)  Stability  bifurcation  diagram  Vx  vs.  o(ilIGJ 

O  Stable,  □  LCO,  A  Unstable 


In  the  absence  of  bending  stiffness  uncertainty,  and  over  a  small  range  of  torsion  stiffness 
uncertainty,  the  wing  remains  stable  until  aaj/GJ  =  0.035,  at  which  the  wing  experiences  flutter 
for  airflow  speed  120  m/s,  above  that  level  the  wing  is  unstable.  Figure  5(b)  shows  the  stability 
bifurcation  diagram  V„-  aGJIGJ .  It  is  seen  that  the  stable  region  in  the  presence  of  torsion 
stiffness  uncertainty  is  smaller  than  the  stable  region  in  the  presence  of  bending  stiffness 
uncertainty.  This  demonstrates  the  significant  influence  of  torsion  stiffness  uncertainty  on  the 
stability  of  the  wing. 

Figure  6  shows  the  influence  of  both  bending  and  torsion  uncertainties  on  the  stability  of  the 
wing.  Figure  6  reveals  the  stability  boundaries  for  the  bending-torsion  combination  uncertainties 
(  geax  /EIx-ogj/  GJ )  at  different  flow  speeds. 
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aa/GJ 


Figure  6.  Stability  boundaries  <j£/j  /  EIX  vs.  aGJ  IGJ  at  different  flow  speeds. 


VII.  MONTE  CARLO  SIMULATION 


In  order  to  define  the  limitation  of  the  perturbation  method,  the  Monte  Carlo  simulation  will  be 
performed.  First  consider  the  deterministic  equations  of  motion  for  one  element  as  given  by 
equation  (5).  Uncertainties  will  be  introduced  in  bending  stiffness,  which  is  described  by 
equation  (6a).  Assembling  the  elemental  matrices  given  in  the  equation  (5)  and  introducing  the 
aerodynamic  loads  given  by  expression  (20),  gives  the  system  equations  of  motion: 


(27) 


MU  +  KU  = 


Using  the  coordinate  transformation  given  by  equation  (23),  equation  (27)  takes  the  form 


(28) 


Equation  (28)  will  be  integrated  numerically  once,  based  on  the  assumption  that  the  response  is 
ergodic.  The  Monte  Carlo  simulation  confirms  the  perturbation  results  in  predicting  the  flutter 
boundary  for  low  values  of  bending  stiffness  variance.  Figure  7(a)  shows  a  comparison  of 
perturbation  (shown  by  solid  curves)  and  Monte  Carlo  simulation  (shown  by  dashed  curves) 
results  for  bending  and  torsion  variances  for  flow  speed  120  m/s,  and  stiffness  uncertainty 
parameters  (tFJi  /EIx  =0.1,  a^/GJ  =  0.0.  Figures  7(b)  show  another  set  for  airflow  speed  of  140 

m/s  and  bending  uncertainty  level  crEI/Eix=  0.015  and  zero  torsion  uncertainty.  Figure  7(b) 
reveals  good  agreement  between  perturbation  and  Monte  Carlo  simulation  results. 
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Figure  7.  Time  history  records  of  bending  and  torsion  response  variance  as  estimated  from  perturbation  method, 
- '  and  Monte  Carlo  Simulation, - ,  for  aa=  5°,  oGJ  IGJ  =  0 ,  qbtn  (0)  =  -0.05 ,  qhen  (0)  =  0 , 

hor  (o)  =  0.01 ,  q,or  (0)  =  0 

(a)  Vn  =  1 20  m  /  s ,  aE,jWx-  0.1.  (b)  F^MOm/s,  <seiJETx  =  0.015. 
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Figure  8.  Stability  bifurcation  diagram  Vx  vs.  ct£V  /£/,  for  Monte  Carlo  Simulation 
O  Stable,  □  LCO,  A  Unstable 

The  bifurcation  diagram  shown  in  Figure  8,  obtained  from  Monte  Carlo  simulation,  is  identical 
with  the  one  shown  in  Figure  5  estimated  by  perturbation  method.  This  demonstrates  that  the 
perturbation  method  is  very  accurate  in  the  prediction  of  the  flutter  boundary  and  the 
establishment  of  stability  boundaries  in  the  presence  of  bending  uncertainties  and  absence  of 
torsion  uncertainties.  The  dependence  of  bending  and  torsion  variances  on  bending  uncertainties 
is  shown  in  Figure  9  for  different  values  of  airflow  speed. 


<*ElJElx  <*ElJEIx 


(a) 

Figure  9.  Comparison  of  perturbation  and  Monte  Carlo  simulation 

(a)  bending,  and  (b)  torsion 


(b) 

response  variances: 


-© - Perturbation  Method; 


—  A - Monte  Carlo  Simulation 
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Figure  10.  Time  history  records  of  bending  and  torsion  response  variances  as  estimated  from  perturbation  method, 
— -  Monte  Carlo  Simulation, - ,  for  ota=5°,  o El%  / EIX  =  0 ,  ^„(0)  =  -0.05,  ^n(0)  =  0,^(0)  =  0.01 , 

<7„,(0)  =  0.  (a)  K,  =120m/s,  /GJ  =  0.05  ,  (b)  V„  =  140m/ s,  IGJ  =  0.01 . 


Figures  10(a)  and  10(b)  show  the  Monte  Carlo  simulation  time  history  response  variances  for 
two  values  of  airflow  speed,  120m/s  and  140m/s,  and  two  different  levels  of  torsion  stiffness 
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uncertainty,  ggj/GJ  =  0.05  and  0.01,  respectively.  The  comparison  shown  in  Figure  10(a)  reveals 
a  poor  approximation  of  the  perturbation  method  as  compared  to  Monte  Carlo  simulation  for 
flow  speed  1 20  m/s  and  aE,x  I  Elx=  0 ,  aGJ  /  GJ  =  0.05 . 

Note  that  the  Monte  Carlo  simulation  yields  LCOs  while  the  perturbation  method  gives  unstable 
response.  Better  correlation  between  the  two  methods  is  obtained  at  higher  speeds  but  lower 
uncertainty  levels  as  demonstrated  in  Figure  10(b)  for  flow  speed  140  m/s,  ge,x/eFx  =  0, 

ggj  IGJ  =  0-01  • 


Figure  1 1.  Stability  bifurcation  diagram  Vx  vs.  ggj/GJ  for  Monte  Carlo  Simulation 
O  Stable,  □  LCO,  A  Unstable 


Figure  12.  Perturbation  method  convergence 
- © - Perturbation  Method;  — A - Monte  Carlo  Simulation 

The  stability  bifurcation  diagram  Vm  -  <sGJ  IGJ  obtained  from  Monte  Carlo  simulation  is  shown  in 
Figure  11.  Comparing  the  bifurcation  diagrams  from  Figures  5(b)  and  11,  one  observes  a  good 
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correlation  in  the  prediction  of  flutter  boundary  up  to  flow  speed  130  m/s  and  ogj/GJ  =  0.03.  As 
the  flow  speed  decreases  and  torsion  stiffness  uncertainty  level  increases,  the  perturbation 
method  loses  the  ability  to  predict  the  flutter  boundary.  For  example,  at  flow  speed  120  m/s,  the 
flutter  onset  is  obtained  form  Monte  Carlo  simulation  at  aaj/GJ  =  0.05,  while  from  perturbation 
method  at  oGJ  IGJ  =  0.035 .  Figure  12(a)  and  12(b)  shows  the  dependence  of  bending  and  torsion 
variances  on  the  torsion  stiffness  uncertainty  level  for  different  values  of  airflow  speed  indicated 
on  each  curve.  Both  diagrams  reveal  the  validity  of  the  perturbation  method  for  a  low  level  of 
stiffness  uncertainty. 


VIII.  CONCLUSIONS 

The  influence  of  uncertainties  of  bending  and  torsion  stiffness  parameters  on  the  flutter  behavior 
of  an  aeroelastic  wing  is  examined.  The  uncertainties  are  modeled  using  a  modified  first  order 
stochastic  perturbation  method  together  with  a  truncated  Karhunen-Loeve  expansion  instead  of 
Taylor  series.  However,  the  Taylor  series  is  used  for  displacement  vector  expansion.  The  results 
of  the  perturbation  approach  are  compared  with  those  predicted  by  Monte  Carlo  simulation  and 
the  comparison  revealed  good  correlation  for  low  values  of  stiffness  uncertainty  levels.  As 
uncertainty  level  increases,  the  perturbation  method  loses  the  accuracy.  For  the  prediction  of 
LCO,  the  perturbation  method  is  very  accurate  for  all  levels  of  bending  stiffness  uncertainty 
examined,  but  the  method  loses  its  accuracy  at  upper  levels  of  torsion  stiffness  uncertainty.  The 
stability  boundary  in  the  flow  speed  versus  stiffness  uncertainty  reveals  the  appearance  of  LCOs 
just  below  the  flutter  speed  boundary.  Further  increase  of  uncertainty  level  produces  instability. 
The  uncertainties  in  torsion  stiffness  induce  a  greater  disturbance  in  the  system.  A  smaller  level 
of  torsion  stiffness  uncertainty  induces  instability  in  the  system. 
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Appendix 

The  detailed  structures  of  Fy  of  equation  (17) 
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where  a  prime  denotes  a  derivative  with  respect  to  y  and  Y,  are  the  finite  element  shape 
functions  given  by 
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CHAPTER  2: 


THE  INFLUENCE  OF  JOINT  RELAXATION 
ON  THE  FLUTTER  OF  AEROELASTIC  STRUCTURES: 


Part  1:  Experimental  Investigation 

ABSTRACT 

The  influence  of  boundary  conditions  relaxation  on  two-dimensional  panel  flutter  is  studied  in 
the  presence  of  in-plane  loading.  The  boundary  value  problem  of  the  panel  involves  time- 
dependent  boundary  conditions  that  are  converted  into  autonomous  form  using  a  special 
coordinate  transformation.  Galerkin’s  method  is  used  to  discretize  the  panel  partial  differential 
equation  of  motion  into  six  nonlinear  ordinary  differential  equations.  The  influence  of  boundary 
conditions  relaxation  on  the  panel  modal  frequencies  and  LCO  amplitudes  in  the  time  and 
frequency  domains  is  examined  using  the  windowed  short  time  Fourier  transform  and  wavelet 
transform.  The  results  revealed: 

•  The  relaxation  and  system  nonlinearity  are  found  to  have  opposite  effects  on  the  time 
evolution  of  the  panel  frequency. 

•  Depending  on  the  system  damping  and  dynamic  pressure,  the  panel  frequency  can 
increase  or  decrease  with  time  as  the  boundary  conditions  approach  the  state  of  simple 
supports. 

The  effects  of  various  bolt  preloads,  viscoelasticity  and  external  applied  static  and  dynamic  loads 
on  bolt  load  relaxation  in  a  carbon/epoxy  composite  bolted  joint  have  been  studied.  Both 
phenomenological  modeling  and  finite  element  analysis  (FEA)  of  bolt-connected  three-point 
bending  specimens  were  employed  in  the  studies.  The  experimental  measurements  showed  that: 

•  Relaxation  of  1 .25%  -  4.25%  over  a  period  of  30  hours  was  observed  depending  on  the 
initial  preload  and  applied  external  loads. 

•  It  was  observed  that  for  any  magnitude  of  external  load  the  bolt  load  relaxation  decreases 
with  increasing  initial  preload.  These  findings  emphasize  the  importance  of  the 
magnitude  of  the  preload. 

•  Only  about  1/3  of  the  relaxation  in  composite  specimens  is  due  to  viscoelastic  behavior  of 
the  polymer  matrix  in  the  composite,  and  the  remaining  2/3  of  the  relaxation  is  due  to 
other  mechanisms  such  as  bolt  thread  slip,  plasticity  and/or  external  excitation. 

State-of-the  Art 

The  purpose  of  the  present  investigation  is  to  conduct  experiments  and  develop 
phenomenological  models  of  relaxation  in  composite  joints  under  external  static  and  dynamic 
loads.  An  attempt  was  also  made  to  indirectly  estimate  the  effect  of  viscoelasticity  on  relaxation 
by  comparing  experimental  relaxation  curves  for  a  composite  bolted  joint  with  those  of  a  steel 
joint,  which  does  not  exhibit  viscoelastic  behavior  at  room  temperature.  Phenomenological 
models  were  derived  from  the  experimental  results  based  on  regression  analysis.  Three- 
dimensional  FEA  models  for  the  composite  bolted  joints  including  the  time-dependent  material 
properties  under  3  point  bending  were  developed  and  a  quasi-elastic  approach  was  used  to 
estimate  the  bolt  load  relaxation  for  static  and  dynamic  loads  and  for  various  values  of  preloads. 
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The  relevant  previous  work  reported  in  the  literature  addresses  the  design,  analysis  and 
relaxation  of  bolted  joints,  and  is  reviewed  in  the  next  two  subsections. 

Design,  Analysis  and  Characterization  of  Joints 

The  design  of  structural  systems  involves  elements  that  are  connected  through  adhesive  bonds, 
bolts,  rivets,  pins,  and  various  combinations  of  these  connections.  Joints  and  fasteners  are  used  to 
transfer  loads  from  one  structural  element  to  another.  In  composite  structures,  three  types  of 
joints  are  commonly  used,  namely,  mechanically  fastened  joints,  adhesively  bonded  joints,  and 
hybrid  mechanically  fastened/adhesively  bonded  joints.  General  aspects  of  design,  analysis  and 
characterization  of  composite  joints  have  been  discussed  by  Camanho  and  Mathews  [1],  Erki  [2], 
Jones  [3],  Wang,  et  al.  [4],  Rastogi,  et  al.  [5,6,7],  and  Hart-Smith  [8-9].  Adhesive  joints  depend 
on  the  size  of  the  parts  to  be  joined  and  the  amount  of  overlap  required  to  carry  the  load. 
Adhesive  joints  are  often  acceptable  for  secondary  structures,  but  are  generally  avoided  in 
primary  structures  on  account  of  strength,  temperature  and  moisture  effects,  and  reliability. 
Bolted  joints  or  hybrid  bolted/bonded  joints  are  still  the  dominant  fastening  mechanisms  used  in 
joining  of  primary  structural  parts  for  advanced  composites. 

The  complex  behavior  of  connecting  elements  plays  an  important  role  in  the  overall  dynamic 
characteristics  such  as  natural  frequencies,  mode  shapes,  and  response  characteristics  to  external 
excitations.  The  joint  represents  a  discontinuity  in  the  structure  and  results  in  high  stresses  that 
often  initiate  joint  failure.  Ibrahim  and  Pettit  [10]  presented  an  overview  of  the  problems 
pertaining  to  structural  dynamics  with  bolted  joints.  The  contact  forces  are  not  ideally  plane  due 
to  manufacturing  tolerances.  Furthermore,  the  initial  forces  will  be  redistributed  non-uniformly 
in  the  presence  of  lateral  loads.  Under  environmental  dynamic  loading,  the  preload  in  joints 
experiences  some  relaxation  that  results  in  time  variation  of  the  structure’s  dynamic  properties. 
Most  of  the  reported  studies  have  focused  on  the  energy  dissipation  of  bolted  joints,  linear  and 
nonlinear  identification  of  the  dynamic  properties  of  the  joints,  parameter  uncertainties  and 
relaxation,  and  active  control  of  the  joint  preload. 

Although  utilized  quite  extensively,  bolted  joints  in  laminated  composites  are  not  well 
understood.  Packman,  et  al.  [11],  Schulz,  et  al.  [12]  and  Lin  and  Jen  [13]  presented  experimental 
and  statistically-based  investigations  of  ultimate  strengths  of  bolted  joints  for  laminated 
composites.  Stress  concentrations  at  hole-edges  in  advanced  composites  can  be  very  high  (Jones, 
[3])  and  joint  efficiencies  are  often  as  low  as  50%  (Hart-Smith,  [14],  The  design  of  composite 
fasteners  subjected  to  transverse  loads  was  considered  by  Running,  et  al.  [15].  Rosner  and 
Rizkalla  [16,17]  experimentally  and  analytically  examined  the  behavior  of  bolted  connections  in 
composite  materials  used  for  civil  engineering  applications.  A  design  procedure  was  introduced 
to  account  for  material  orthotropy,  pseudo-yielding  capability,  and  other  factors  that  influence 
bolted-connection  behavior. 

Andreasson,  et  al.  [18],  Menendez  and  Guemes  [19],  and  Li,  et  al.  [20]  analytically  and 
experimentally  studied  the  tensile  response  and  failure  of  bolted  or  riveted  joints  made  from 
carbon  fiber-reinforced  plastic  (CFRP).  They  observed  that  the  dynamic  behavior  of  composite 
joints  is  much  more  complicated  than  its  behavior  for  the  quasi-static  condition  due  to  strain  rate 
and  inertia  effects.  The  effect  of  clamping  on  the  tensile  strength  of  composite  plates  with  a  bolt- 
filled  hole  was  studied  by  Horn  [21],  Hung  [22],  Hung,  et  al.  [23,24],  and  Yan,  et  al.  [25].  Early 
work  by  Van  Siclen  [26]  on  composite  bolted  joints  made  of  graphite/epoxy  involved  the 
development  of  joint  allowables  and  alternating  reinforcing  concepts  for  increasing  joint 
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strength.  Semi-empirical  procedures  combined  with  experimental  data  were  used  to  predict  the 
joint  strength  as  a  function  of  basic  laminate  properties  and  joint  geometrical  parameters.  Studies 
reported  in  the  literature  address  the  importance  of  joint  allowables  [27,  28],  selection  of 
optimum  ply  stacking  sequence  [8],  washer  size,  effects  of  clamping  force  or  bolt  preload  on 
joint  strength  and  failure  modes  [9,29,30].  Yet,  some  of  these  questions  remain  unanswered.  For 
example,  how  much  bolt  preload  should  be  applied  to  composite  joints  to  minimize  relaxation 
and  what  other  factors  have  the  most  influence  on  relaxation  in  such  joints? 

Some  designers  use  the  existing  design  code  MSFC-STD-486B  from  the  NASA  Marshall  Space 
Flight  Center  [31],  which  was  developed  originally  for  metallic  joints.  This  code  recommends 
that,  for  typical  preloaded  structural  assemblies  in  tension,  the  bolt  preload  should  be  less  than  65 
percent  of  the  tensile  yield  strength  of  the  fastener,  while  the  corresponding  bolt  preload  for  lap- 
shear  joints  should  be  less  than  60  percent  of  the  fastener  yield  strength.  However,  the  above 
preload  values  are  directly  related  to  the  tension  in  the  fastener  but  not  the  structural  members  to 
be  connected.  In  the  case  of  composite  joints,  only  half  of  the  bolt  preload  value  recommended 
for  steel  joints  has  been  suggested  for  any  given  bolt  size  [32].  Thomas  and  Zhao  [33]  tested 
composite  joints  made  of  graphite/epoxy  with  different  thicknesses  and  bolt  diameters  and  found 
that  preload  limits  as  specified  by  MSFC-STD-486B  are  acceptable,  but  recommended  that  a 
torque-tension  test  be  performed  while  simultaneously  monitoring  the  energy  levels  using 
acoustic  emissions  to  prevent  damage  to  the  joint.  The  reason  for  this  is  that  mechanical 
fasteners  such  as  bolts  and  rivets  typically  generate  clamping  forces  in  the  through-the-thickness 
(TTT)  direction  of  composite  structures,  perpendicular  to  the  direction  of  fiber  reinforcement. 
Such  structures  are  known  to  have  poor  properties  in  the  TTT  direction  and  are  susceptible  to 
damage  and  failure.  This  is  particularly  important  for  viscoelastic  behavior,  because  the  time- 
dependent  viscoelastic  behavior  of  polymer  matrix  composites  is  typically  most  evident  in  the 
matrix-dominated  TTT  direction. 

Cooper  and  Turvey  [34]  investigated  the  effects  of  bolt  tightening  torque  on  failure  loads,  critical 
joint  geometry  and  joint  stiffness  in  single  bolt  double  lap  tension  joints  in  pultruded  glass 
reinforced  plastics  (GRP)  sheet  material.  It  was  observed  that  as  the  bolt  tightening  torque  was 
increased  from  3  N-m  (lightly  torqued)  to  30  N-m  (heavily  torqued)  the  bearing  failure  load 
increased  by  about  50%  and  the  critical  end  distance,  determined  by  the  E/D  ratio  (where  E  is  the 
edge  distance,  and  D  is  the  bolt  diameter)  which  ensures  bearing  failure,  was  also  found  to 
increase,  but  the  stiffness  of  the  joint  remained  unaffected.  It  was  also  shown  by  Lim,  et  al.  [35] 
that  for  a  glass/epoxy  bolted  composite  joint  with  a  stacking  sequence  of  [±9/0g]s  static  shear 
strengths  increased  by  about  30%  when  clamping  pressure  was  varied  from  0  to  60  MPa.  In 
addition,  the  specimens  with  lower  preloads  failed  at  80,000  cycles  whereas  specimens  with 
higher  preloads  lasted  up  to  3  million  cycles.  In  all  of  the  above  studies  the  clamp-up  forces  were 
selected  randomly,  the  effects  of  external  loads  on  relaxation  were  not  considered  and  in  some 
cases  the  relaxation  in  the  bolted  joint  was  not  monitored. 

Bolt  Preload  Relaxation  Under  External  Excitation 

There  are  several  factors  affecting  relaxation  in  bolted  metallic  joints,  which  are  well 
documented  by  Bickford  [36].  A  fastener  subjected  to  vibration  will  not  lose  all  pre-loads 
immediately.  First  there  will  be  a  slow  loss  of  pre-load  caused  by  some  of  the  relaxation 
mechanisms.  Vibration  will  increase  relaxation  because  wear  and  hammering  will  take  place 
during  vibration.  After  sufficient  pre-load  is  lost,  friction  forces  drop  below  a  critical  level  and 
the  nut  actually  starts  to  back  off  and  shake  loose.  In  this  case,  the  joints  will  not  resemble  the 
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ideal  boundary  conditions  but  will  involve  uncertainties.  With  higher  initial  pre-load,  longer  or 
more  severe  vibration  is  required  to  reduce  pre-load  to  the  critical  point  at  which  back-off  occurs. 
In  fact,  in  some  circumstances,  if  pre-load  is  high  enough  to  start  with;  nut  back  off  will  never 
take  place.  It  is  not  known  if  these  observations  apply  to  composite  bolted  joints.  Jiang  et  al.  [37] 
conducted  experimental  investigations  and  finite  element  analysis  (FEA)  on  relaxation  in  steel 
joints  under  cyclic  loading.  Depending  on  the  loading  conditions  and  preload,  the  loss  of 
clamping  force  ranged  from  a  few  percent  to  40%  after  200  cycles.  Three-dimensional  FEA  with 
an  advanced  cyclic-plasticity  model  was  implemented  to  predict  the  clamping  force  reduction. 

Schmitt  and  Horn  [38],  and  Horn  and  Schmitt  [39]  studied  the  relaxation  process  in  bolted 
thermoplastic  composite  joints.  From  single  lap-shear  tests  of  bolted  thermoplastic  composite  lap 
joints  involving  viscoelastic  effects,  Horn  and  Schmitt  [39]  observed  that  by  increasing  clamp-up 
force  the  joint  bearing  failure  load  increased  by  as  much  as  28%.  Clamp-up  force  relaxation  was 
also  monitored  and  6%  -  16%  relaxation  from  the  initial  clamp-up  force  was  observed  for  short 
duration  (1400  hours)  and  long  duration  (100,000  hours),  respectively.  The  rate  of  relaxation  was 
predicted  to  be  as  high  as  37%  at  an  elevated  temperature  of  250°  F  for  the  duration  of  100,000 
hours.  Zhao  and  Gibson  [40]  showed  that  compressive  clamping  stress  relaxed  by  1 8%  and  1 5% 
in  E-glass/epoxy  beams  with  and  without  polymeric  interleaves,  respectively,  for  a  period  of  50 
hours,  whereas  the  corresponding  relaxation  in  an  aluminum  beam  was  negligible. 

Using  finite  element  analysis  (FEA),  Shivakumar  and  Crews  [41]  predicted  bolt  force  relaxation 
of  up  to  31%  at  room  temperature  for  a  20  year  duration  in  a  double  lap  joint  made  of 
graphite/epoxy.  It  was  also  observed  that  the  rate  of  relaxation  increases  with  increasing 
temperature  and  moisture  levels.  But  these  findings  were  for  a  randomly  selected  clamp-up  force 
and  neither  the  effect  of  different  clamp-up  forces  nor  the  effects  of  external  loads  on  relaxation 
were  considered. 


EXPERIMENTAL  INVESTIGATION 


Experimental  Setup 

Figure  1  shows  the  configuration  of  the  specimen  in  3-point  bending.  The  design  of  the  specimen 
support  was  selected  in  such  a  way  to  maintain  combined  bending  and  shear  loading.  The  joints 
were  made  from  unidirectional  carbon/epoxy  prepreg  tapes  (Prepreg  type:  P2254-20-305  T800 
carbon/epoxy)  from  Toray  Composites  America,  Inc.  A  total  of  48  layers  of  prepreg  were  laid  on 
top  of  each  other  in  the  mold  along  with  bleeder  cloth  and  release  fabric.  The  mold  assembly  was 
cured  in  the  molding  chamber  of  a  TMP  autoclave-style  vacuum  press,  and  the  cured  thickness 
of  the  sample  plate  was  found  to  be  7.75  mm.  Two  rectangular  beams  were  cut  from  the  plate 
and  machined  to  the  required  dimensions  as  shown  in  Figure  1.  Joint  clamping  force  was 
generated  by  an  instrumented  hexagonal  head  steel  bolt  with  built-in  strain-gages  (Strainsert® 
Model  SXS-FB,  9.525  mm  or  3/8  inch  in  diameter).  The  driving  force  for  stress  relaxation  is  the 
imposed  strain,  as  governed  by  the  initial  strained  length  of  the  bolt.  This  is  a  relaxation  problem, 
not  a  creep  problem.  The  creep  strain  under  applied  stress  depends  on  the  stress  level,  but  the 
stress  relaxation  depends  on  the  imposed  strain,  which  should  be  independent  of  washer  size. 
Therefore,  only  one  washer  size  (19  mm  diameter)  was  used  to  clamp  the  composite  specimen 
from  either  side. 

Two  separate  data  acquisition  systems  were  used,  one  for  the  instrumented  bolt  and  another  to 
control  the  Enduratec  servo-pneumatic  testing  machine.  The  beam  load  and  displacement  values 


were  monitored  by  transducers  which  were  built  in  to  the  Enduratec  machine.  The  bolt  preload 
selection  was  based  on  the  maximum  tensile  load  recommended  by  the  manufacturer 
(Strainsert®).  Bolt  preloads  selected  for  this  study  were  4200  N,  5050  N,  6700  N  and  7850  N, 
which  correspond  to  12.5%,  15%,  20%  and  23.5%,  respectively,  of  the  maximum  manufacturer- 
recommended  tensile  load  in  the  instrumented  bolt.  The  corresponding  through-the-thickness 
compressive  stresses  generated  in  the  composite  specimen  were  approximately  25%,  29%,  38% 
and  45%,  respectively,  of  the  estimated  transverse  compressive  strength  of  the  composite.  After 
applying  an  initial  preload  to  the  composite  bolted  joint,  the  bolt  load  was  monitored  for  a  period 
of  30  hours  for  the  following  cases:  1)  bolt  preload  in  the  absence  of  external  beam  load,  2)  bolt 
preload  plus  a  static  applied  load  of  250  N  at  beam  midspan  (combination  of  ramp  and  dwell), 
and  3)  bolt  preload  plus  a  dynamic  load  of  250  N  amplitude  and  frequencies  of  1  Hz  and  5  Hz  at 
beam  midspan  as  shown  in  Figure  2(a).  Only  one  composite  joint  specimen  was  manufactured 
and  tested.  However,  some  of  the  relaxation  tests  were  repeated  with  the  same  specimen  under 
different  loading  conditions  (for  example:  5050  N  bolt  preload  at  5  Hz  and  7850  N  preload  only 
condition).  It  was  found  that  the  relaxation  was  reproducible  and  the  differences  in  the  relaxation 
from  one  test  to  another  were  less  than  0.5%.  A  stability  check  on  the  data  acquisition  system 
was  carried  out  by  monitoring  the  signal  from  the  instrumented  bolt  for  a  period  of  30  hours 
without  the  application  of  load  to  the  bolts.  The  drift  in  the  signal  for  the  above  period  was  less 
than  1%.  Also,  the  drift  in  the  signal  was  taken  into  account  when  the  percentage  relaxation  was 
calculated. 

Static  and  dynamic  load  experiments  were  carried  out  on  a  3-point  bend  fixture  mounted  on  the 
Enduratec  servo-pneumatic  testing  machine  as  shown  in  Figure  2(b).  For  the  first  case, 
experiments  were  also  conducted  on  a  steel  joint  having  the  same  dimensions  as  those  of  the 
composite  joint.  Since  the  relaxation  in  the  steel  joint  was  expected  to  be  mainly  due  to  thread 
slip  and/or  other  mechanisms  rather  than  viscoelastic  material  behavior,  these  tests  were  done  in 
order  to  better  interpret  the  results  of  the  composite  joint  tests.  The  bolt  load  was  monitored  for  a 
period  of  30  hours  using  the  same  3-point  bend  set-up  and  the  instrumented  bolt  for  preloads 
mentioned  above.  All  the  above  experiments  were  conducted  without  any  thread  seal.  In  an 
attempt  to  gain  some  insight  into  the  importance  of  slip  in  the  bolt  threads  as  one  of  the  possible 
relaxation  mechanisms,  another  set  of  experiments  was  conducted  in  the  absence  of  any  external 
applied  load  with  Teflon  tape  thread  seal  wrapped  around  the  threads.  The  Teflon  thread  seal 
(Polytetrafluoroethylene,  or  PTFE)  meets  standard  MIL  T-27730A.  The  bolt  load  was  monitored 
for  the  duration  of  30  hours. 

EXPERIMENTAL  RESULTS  AND  PHENOMENOLOGICAL  MODELING 
Effect  of preload  on  relaxation  in  composite  and  steel  joints 

Experimental  measurements  of  bolt  load  relaxation  for  composite  and  steel  bolted  joints  with 
initial  preloads  of  4200  N,  5050  N,  6700  N  and  7850  N  and  no  applied  beam  load  are  shown  in 
Figure  3.  This  figure  represents  the  time  evolution  of  the  preload  (represented  by  the  normalized 
load  parameter,  P  =  P(t)/P0 ,  where  P0  is  the  initial  bolt  preload  and  P(t)  is  the  relaxed  preload 

at  time  t .  It  is  observed  that  for  the  composite  bolted  joint,  the  magnitude  of  the  relaxation 
decreases  with  increasing  bolt  preload.  For  steel  joints,  the  magnitude  of  the  relaxation  is 
approximately  the  same  for  the  preloads  selected  and  the  relaxation  for  a  period  of  27  hours  is 
less  than  that  of  the  composite  joint.  Relaxation  of  approximately  4%  in  the  bolt  preload  (4200 
N)  for  the  composite  joint  was  observed  for  a  period  of  30  hours,  but  the  relaxation  decreased  to 
approximately  3  %  for  higher  preloads  (7850  N).  Figure  4  reveals  that  the  preload  relaxation 
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magnitude  in  composite  joints  increases  when  bolted  with  a  thread  seal,  indicating  that  a 
significant  portion  of  the  relaxation  observed  is  due  to  slip  in  the  threads  and/or  viscoelastic 
relaxation  in  the  Teflon  thread  seal.  However,  it  is  again  observed  that  the  bolt  load  relaxation 
decreases  with  increasing  initial  bolt  preload. 

Effect  of  applied  static  loading 

Figure  5  shows  the  effect  of  static  beam  load  on  bolt  load  relaxation  in  the  composite  bolted  joint 
under  the  3  point  bend  setup.  It  is  seen  that  the  magnitude  of  the  relaxation  again  decreases  with 
increasing  bolt  preload.  However,  when  compared  with  the  preload  only  condition,  the  trend  of 
decreasing  bolt  load  relaxation  with  increasing  bolt  preload  only  becomes  clear  for  higher 
preloads  (6700  N  and  7850  N)  and  the  result  for  the  lowest  preload  (4200  N)  seems  inconsistent 
with  that  trend. 

Effect  of  upplied  dynamic  loading 

Figures  6  and  7  show  the  preload  relaxation  plots  under  dynamic  applied  loading  on  the 
beam  F(t)  =  F0  +  AFsin  2 nft ,  where  F0  is  the  static  component,  A F  =  250  N  is  the  amplitude  of 
the  sinusoidal  component,  and  /  is  the  frequency  taken  as  1  Hz  or  5  Hz.  It  is  seen  that  under  the 
excitation  frequency  of  1  Hz,  the  magnitude  of  preload  relaxation  decreases  with  increasing  bolt 
preload  and  is  nearly  the  same  for  higher  preloads  (6700  N  and  7850  N).  However,  as  the 
frequency  of  the  external  beam  load  is  increased  to  5  Hz,  the  magnitude  of  the  relaxation  tends  to 
increase  as  shown  in  Figure  7.  For  lower  preloads  (5050  N)  the  magnitude  of  relaxation  in  the 
bolt  preload  exceeds  that  of  the  higher  preload  cases  during  the  initial  time  interval  of  the 
experiment,  however,  this  trend  reverses  for  longer  times.  The  increase  in  the  magnitude  of 
relaxation  at  higher  frequencies  of  excitation  and  longer  times  may  be  due  to  the  increase  in  the 
number  of  cycles  as  reported  in  Bickford  [36].  This  is  associated  with  a  corresponding  increase 
in  frictional  heating.  The  temperature  at  the  interface  of  the  lap  joint  was  measured  before  and 
immediately  after  the  experiment  with  a  general  purpose  type  K  (Chromel  /  Alumel) 

thermocouple  and  was  found  to  increase  by  3 °C,  and  1°  C,  respectively,  when  the  composite 
bolted  joint  was  preloaded  to  6700  N  and  excited  with  5  Hz  and  1  Hz  frequency,  respectively,  for 
the  duration  of  30  hours.  When  the  relaxation  in  bolt  load  due  to  applied  dynamic  load  at  5  Hz  is 
compared  with  those  due  to  applied  static  loads,  it  is  found  that  for  lower  preloads  (5050  N)  the 
applied  static  load  increases  the  bolt  preload  relaxation  whereas  for  higher  bolt  preloads  (6700  N 
and  7800  N)  the  bolt  load  relaxation  decreases. 

PHENOMENOLOGICAL  MODELING 

The  experimental  bolt  force  relaxation  data  in  Figures  3-7  can  be  adequately  described  by  an 
inverse  power  law  of  the  type  suggested  by  Shivakumar  and  Crews  [41].  In  normalized 
dimensionless  form,  this  equation  can  be  written  as 

P(t)  _  1 

P0  1  +  Kt"  (1) 

where 

P0  =  initial  bolt  preload 

P(t)=  relaxed  bolt  preload  at  time  / 

K  =  constant  in  (hours)  'N 
N  =  dimensionless  exponent 
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t  =  time  in  hours 

For  the  data  in  Figures  3-7,  the  constants  K  and  N  in  Equation  (1)  were  evaluated  using  the  least- 
squares  regression  analysis  in  MATLAB  [42]  at  the  99%  confidence  level  and  are  tabulated  in 
Tables  (la  -  If)  for  the  different  preloads  and  beam  loads.  For  each  case  the  coefficient  of 
multiple  determination,  R2,  is  also  given.  This  coefficient  is  a  statistical  measure  of  the  goodness 
of  the  fit  of  the  equation  to  the  data,  where  R2  =  1.0  corresponds  to  a  perfect  fit.  It  should  be 
noted  that  the  constants  K  and  N  listed  in  Table  (la  -  If)  are  valid  only  up  to  30  hours,  which  is 
the  maximum  test  duration.  It  is  seen  that  there  are  no  consistent  trends  in  the  variations  of  AT  and 
N  with  increasing  preload,  but  the  combined  effect  of  the  variations  in  these  two  quantities  is 
such  that  increasing  preload  reduces  the  relaxation. 

VISCOELASTIC  ANALYSIS  AND  MATERIAL  CHARACTERIZATION 

The  three-dimensional  (3-D)  elastic  lamina  properties  needed  for  the  finite  element  analysis  were 
calculated  from  composite  micromechanics  equations  [43]  (for  example:  the  rule  of  mixtures  to 
calculate  the  longitudinal  modulus  and  Tsai-Hahn  equations  to  calculate  the  transverse  modulus 
and  shear  modulus)  using  the  fiber  and  matrix  properties  together  with  the  following 
assumptions: 

1 .  The  fibers  are  linear  elastic. 

2.  The  matrix  is  linear  viscoelastic,  with  its  creep  compliance  described  by  a  power  law 

3.  The  composite  is  specially  orthotropic  and  transversely  isotropic 

4.  The  viscoelastic  response  depends  only  on  the  time  elapsed  since  application  of  the  load 
(i.e.  the  material  is  assumed  to  be  non  aging). 

The  analysis  consisted  of  three  parts.  First,  since  the  viscoelastic  creep  and  relaxation  data  for  the 
epoxy  resin  used  in  the  prepreg  is  not  available,  Beckwith’s  [44]  measured  linear  viscoelastic 
properties  for  Shell  58-68  epoxy  at  75°  F  were  assumed.  Next,  these  properties  were  used  in  the 
FEA  to  predict  the  bolt  load  relaxation  for  an  epoxy  beam  under  different  loading  conditions 
using  the  quasi-elastic  approach.  Since  ABAQUS  viscoelastic  modeling  capability  is  limited  to 
isotropic  materials,  it  was  necessary  to  validate  the  quasi-elastic  approach  (explained  in  the  next 
section).  Beckwith’s  [44]  creep  test  results  were  extrapolated  out  to  50  hrs  from  the  available 
data,  by  using  the  empirical  power  law  equation  for  creep  compliance: 

D{t)  =  D0  +  ZV”  (2) 


Where,  from  [44], 

D(t)  =  time  dependent  isotropic  creep  compliance  of  matrix 
D0  =  initial  elastic  compliance  of  matrix  =  2.726  x  10'4  (MPa)'1 
D/~  creep  coefficient  of  matrix  =  1 .0  x  1  O'3  'MPa)'1  (min)'" 
t  =  time  in  minutes 

n  =  dimensionless  creep  exponent  =  0.19 


The  time-dependent  viscoelastic  properties  of  the  composite  joint  were  assumed  to  depend  only 
on  the  time-dependent  properties  of  the  epoxy  matrix  material.  Based  on  the  linear  viscoelastic 
assumption,  a  time-dependent  matrix  modulus,  Em(t),  was  estimated  from  the  following  equation: 


3.(0* 


1 

D(t) 
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(3) 


It  was  assumed  that  the  power  law  constants;  Do,  D /  and  n  at  room  temperature  are  the  same  as 
those  measured  by  Beckwith  for  Shell  58-68  epoxy  [44],  By  using  these  constants  and  the  power 
law,  the  time-dependent  creep  compliance,  D(t),  for  the  epoxy  material  was  calculated.  The  time- 
dependent  Young’s  modulus  of  the  matrix  material,  Em(t),  (see  Table  2)  was  found  from 
Equation  (3).  The  tensile  modulus  data  for  T800  carbon  fibers  was  obtained  from  the  fiber 
manufacturer,  Toray  Composites  America,  Inc.  as  £/  =  294  GPa  and  was  assumed  to  be 
independent  of  time.  Viscoelastic  properties  of  the  lamina  (i.e.,  longitudinal  modulus  Ei(t), 
transverse  moduli  E2(t)  and  E}(t),  Poisson’s  ratios  vi2(t),  v/3(t)  and  »23(t),  longitudinal  shear 
modulus  Gn(t),  and  transverse  shear  moduli  G/3(t)  and  G23(t)  were  estimated  using  elastic  fiber 
properties  and  time-dependent  viscoelastic  resin  properties  through  an  application  of  the  Elastic- 
Viscoelastic  Correspondence  Principle  to  the  micromechanics  equations  [43]  with  a  fiber  volume 
fraction,  Vf  =  0.45  (see  Table  2).  Fiber  volume  fraction  for  the  composite  laminate  was  indirectly 
estimated  using  the  combined  experimental/numerical  technique.  Measured  load-displacement 
response  from  static  3  point  bend  tests  was  compared  with  the  predicted  response  (FE  models) 
initially  assuming  the  fiber  volume  fraction  to  be  0.6.  The  difference  between  the  predicted  and 
measured  response  was  then  minimized  using  fiber  volume  fraction  as  curve  fitting  parameter. 
These  properties  were  then  used  in  the  FEA  to  calculate  the  bolt  load  relaxation. 

Due  to  the  limitations  of  ABAQUS  [45],  the  quasi-elastic  approach  [43]  was  used  to  predict  the 
bolt  load  relaxation  for  the  orthotropic  composite  beams.  In  this  approach,  the  viscoelastic 
solutions  were  approximated  by  a  series  of  elastic  solutions  corresponding  to  different  elastic 
properties  at  different  times,  while  the  stresses  were  assumed  to  be  constant  within  each  time 
increment. 


FINITE  ELEMENT  ANALYSIS 


Numerical  Modeling 

The  purpose  of  the  FEA  is  to  develop  predictive  numerical  models,  and  to  promote  a  more 
meaningful  interpretation  of  the  experimental  results.  Hypermesh®  5.0  [46]  was  used  for  model 
development  and  post-processing,  and  ABAQUS®  V  6.3  Standard  3-D  [45]  was  used  predict  the 
bolt  load  relaxation  using  the  quasi-elastic  analysis.  In  order  to  benchmark  the  quasi-elastic 
prediction  against  the  theoretically  “exact”  viscoelastic  solution  in  ABAQUS  (which  is  only 
applicable  to  isotropic  materials),  comparisons  were  made  between  ABAQUS  viscoelastic 
solution  and  the  quasi-elastic  relaxation  predictions  for  an  un-reinforced  isotropic  epoxy  beam 
[47].  The  only  reason  for  this  part  of  the  study  was  to  compare  the  relaxation  predictions  made 
from  viscoelastic  (exact)  analysis  with  those  from  the  quasi-elastic  analysis.  In  these  studies,  the 
bolt  preload  relaxation  was  predicted  for  different  static  and  dynamic  loads  using  both 
viscoelastic  and  quasi-elastic  approaches.  Relaxation  of  about  4.75%  in  the  bolt  preload  was 
predicted  for  a  period  of  50  hours,  and  the  difference  between  viscoelastic  and  quasi-elastic 
predictions  was  shown  to  be  less  than  2%  [47]. 

With  the  confidence  gained  from  the  benchmarking  study,  a  FEA  model  (Figure  8)  for  the 
composite  bolted  joint  was  developed  using  the  quasi-elastic  approach  only.  A  global-local 
submodeling  technique  was  used  to  model  the  bolted  composite  joint,  where  the  displacements 
around  the  bolted  joint  section  in  the  global  model  (a  one-piece  beam  with  no  bolt)  are  used  to 
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drive  the  refined  local  model  (Figure  8).  This  technique  has  the  advantage  that  more  detailed 
results  in  the  vicinity  of  the  bolted  joint  can  be  obtained  with  fewer  elements  than  with  the  full 
model.  All  models  were  developed  using  type  C3D8  and  C3D6  3-D  solid  elements.  Since  the 
model  is  symmetric  about  the  vertical  midplane  of  the  beam,  only  a  half  model  was  used,  and 
symmetric  boundary  conditions  were  applied.  A  mesh  sensitivity  study  was  performed  for  both 
global  and  local  model  to  ensure  that  the  finite  element  meshes  were  fine  enough  to  give 
satisfactory  results.  For  example,  the  bolt  preload  predicted  by  the  local  model  with  4664 
elements  was  compared  with  the  local  model  with  a  coarser  mesh  having  only  2857  elements, 
and  the  difference  was  found  be  5%.  By  increasing  the  number  of  elements  in  the  local  model  to 
5632  elements  the  model  was  found  to  be  converging  and  the  difference  in  the  predicted  bolt 
load  was  found  be  less  than  1%.  Therefore,  all  of  the  analyses  were  carried  out  using  the  local 
model  with  4664  elements. 

The  threads  in  the  bolt  were  neglected  in  the  FEA  models,  and  the  bolt  was  assumed  to  be  a  solid 
cylinder.  Thus,  possible  bolt  load  relaxation  due  to  plastic  deformation  and/or  thread  slip  in  the 
threads  was  not  included  in  the  models.  The  solid  bolt  simulation  requires  that  contact  surfaces 
be  defined  between  all  the  surfaces  that  are  in  contact,  and  these  surfaces  were  modeled  using  the 
contact  pair  approach  in  ABAQUS.  The  contact  pairs  are  defined  from  free  element  faces.  Since 
the  sliding  between  the  surfaces  was  expected  to  be  small,  the  ‘small  sliding’  option  was  used  in 
all  analyses.  Friction  coefficients  were  set  to  0.2  for  all  contact  surfaces,  as  used  by  Ireman  [48], 
The  coefficient  of  friction  (COF)  is  very  difficult  to  control  and  measure  or  predict  as  it  depends 
on  the  surface  conditions  of  the  joining  parts.  Several  authors  [49-52]  have  measured  COFs 
(both  static  and  kinetic)  of  polymer  materials  without  any  reinforcements,  polymeric  based 
composite  materials  in  contact  with  composites  as  well  as  when  in  contact  with  metals.  COFs 
reported  in  the  literature  range  from  as  low  as  0.096  to  as  high  as  0.74.  The  COF  has  been  shown 
to  be  independent  of  the  applied  normal  load  [49,  50].  Herrington,  et  al.  [49]  found  the  static 
friction  coefficient  to  be  in  the  range  of  0.096  -  0.128  for  a  simple  double  shear  arrangement 
(without  bolts)  with  steel  plates  on  the  outer  surface  of  the  composite  specimen.  Schon  [50,  51] 
observed  that  the  friction  coefficient  initially  for  composite-composite  is  about  0.65  but  only 
0.23  when  the  composite  is  in  contact  with  aluminum  [51].  Xiao  et  al.  [52]  have  measured  the 
COF  to  be  in  the  range  of  0.125  -  0.3,  when  a  steel  pin  is  sliding  against  the  edge  of  a 
unidirectional  carbon  fiber/epoxy  (T800H/3631)  laminate.  It  should  also  be  noted  that  several 
studies  found  in  the  literature  consider  frictionless  pins  or  bolted  models  for  their  analysis.  Since 
“small  sliding”  contact  is  assumed  for  the  current  study,  the  use  of  COF  equal  to  0.2  seems 
reasonable. 

The  beam  load  was  applied  as  a  concentrated  nodal  force  at  midspan.  Pretension  in  the  bolt  was 
applied  in  a  separate  loading  step  by  defining  a  pre-tension  section  in  the  bolts.  Assembly  loads 
were  transmitted  across  the  pre-tension  section  by  means  of  a  pre-tension  node.  Pre-load  was 
applied  by  giving  an  initial  displacement  (in  the  direction  parallel  to  the  bolt  axis)  to  this  node. 
Bolt  preload  was  maintained  by  using  the  fixed  option  under  boundary  conditions,  and  was 
monitored  by  checking  the  total  force  output  on  that  node.  The  model  geometry  and  boundary 
conditions  were  the  same  as  for  the  experiments,  as  shown  in  Figure  1.  For  the  composite  beams, 
both  experiments  and  FEA  (using  the  global-local  model  and  quasi-elastic  approach)  were 
conducted,  but  only  experiments  were  conducted  for  steel  joints. 

Bolt  load  relaxation  was  predicted  for  a  period  of  30  hours  in  composite  bolted  joints  using  a 
quasi-elastic  analysis  and  the  material  properties  listed  in  Table  2  for  a  preload  of  4200  N  under 
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the  following  types  of  loading:  1)  bolt  preloading  in  the  absence  of  external  beam  loading,  2) 
bolt  preloading  in  the  presence  of  a  static  beam  load  of  250  N,  and  3)  bolt  preloading  in  the 
presence  of  dynamic  loading  of  amplitude  250  N  at  a  frequency  of  1  Hz.  The  dynamic  beam  load 
was  applied  using  the  periodic  loading  option  in  ABAQUS  for  5  cycles  under  each  quasi-elastic 
step  and  mean  bolt  load  data  was  recorded  in  each  case 

Numerical  Results 

Figure  9  shows  the  predicted  bolt  load  relaxation  for  composite  beams  under  static  and  dynamic 
beam  loads  using  the  quasi-elastic  approach  with  the  material  properties  listed  in  Table  2.  It  is 
observed  that  the  predicted  bolt  loads  are  shifted  slightly  with  the  application  of  the  static  and 
dynamic  beam  loads,  but  otherwise  the  relaxation  curves  are  not  significantly  affected  by  the 
beam  loads.  Again,  the  mean  bolt  load  is  plotted  for  the  dynamic  analysis. 

Even  though  the  FEA  results  captured  the  bolt  load  relaxation  in  composite  joints  with  and 
without  external  loads.  Figure  10  reveals  that,  in  the  absence  of  external  loading,  there  is  no 
change  in  the  predicted  magnitude  of  bolt  load  relaxation  with  increasing  bolt  preloads.  This 
contradicts  the  experimental  observation  that  relaxation  decreases  with  increasing  bolt  preload. 
There  are  several  possible  reasons  for  this  disagreement. 

First,  the  material  model  does  not  capture  the  viscoelastic  effects  in  the  polymer  matrix  material 
at  the  micromechanical  level.  The  micromechanical  analysis  referred  to  in  the  discussion 
following  Equation  3  was  based  on  “mechanics  of  materials”  type  models  which  do  not  take  into 
account  the  details  of  the  in-situ  stress  and  strain  distributions  in  the  viscoelastic  polymer  matrix. 
A  3-D  finite  element  micromechanics  model  which  takes  into  account  the  micromechanical 
geometry  is  needed  to  accurately  simulate  the  effects  of  such  parameters  as  boundary  conditions 
at  the  bolt-composite  interface  on  the  viscoelastic  relaxation  of  the  composite.  In  the  absence  of 
such  a  detailed  3-D  micromechanics  analysis,  an  analysis  of  the  effects  of  preload  on  the 
macromechanical  volume-averaged  von  Mises  stress  was  conducted.  One  measure  of  the 
principal  stress  differences  and  corresponding  shear  stresses  is  the  von  Mises  stress  given  by 


(4) 


Where,  crp  <x2  and  are  principal  stresses  in  1,  2  and  3  directions,  respectively.  Figure  11 

shows  a  comparison  of  macromechanical  volume-averaged  (or  “effective”)  von  Mises  stresses, 
o\. ,  in  the  immediate  vicinity  of  the  bolt  hole  in  the  composite  material  for  different  preloads. 

One  can  conclude  that  the  contribution  of  the  preload  to  the  total  effective  von  Mises  stresses 
(i.e.  for  preload  +  applied  static  load)  clearly  increases  with  increasing  bolt  preloads.  One  would 
normally  expect  that  the  increase  in  the  von  Mises  stress  would  be  associated  with  a 
corresponding  increase  in  the  tendency  to  exhibit  viscoelastic  behavior.  However,  the 
experimental  observation  indicates  reduced  relaxation  with  increased  bolt  preload.  The  results  in 
Figure  1 1  are  based  on  a  macromechanical  analysis,  and  a  micromechanical  analysis  of  the 
viscoelastic  behavior  of  the  polymer  matrix  material  may  lead  to  a  different  conclusion.  Another 
possible  reason  for  the  FEA  model’s  inability  to  predict  the  effect  of  increasing  bolt  preload  on 
the  relaxation  response  is  that  the  initial  radial  gap  between  the  outside  diameter  of  the  bolt  and 
the  inside  diameter  of  the  bolt  hole  in  the  composite  may  be  closed  as  the  increasing  bolt  preload 
squeezes  the  composite  and  causes  the  hole  size  to  shrink.  This,  in  turn,  could  lead  to  time- 
dependent  boundary  conditions  at  the  inner  surface  of  the  hole  as  the  relaxation  proceeds.  As  the 
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composite  presses  against  the  bolt,  the  effect  would  be  to  increase  the  constraint  on  the 
composite  and  reduce  its  relaxation  response.  Further  work  in  this  area  is  clearly  needed. 

COMPARISONS  AND  VALIDATION  WITH  EXPERIMENTAL  RESULTS 

FEA  results  obtained  for  composite  bolted  joints  with  4200  N  preload  under  preload  only,  with 
static  preload  of  250  N  and  dynamic  beam  load  at  1  Hz  frequency  are  compared  with  the 
experimental  bolt  load  relaxation  in  Figures  3,  5  and  6.  FEA  results  agree  very  well  with  the 
experimental  observations  for  a  Fixed  preload  of  4200  N.  The  model  not  only  captures  the 
reduction  in  the  bolt  load  during  the  first  few  hours,  but  also  the  general  tendency  of  decay  in 
bolt  load  with  increasing  loading  cycles.  However,  as  pointed  out  in  the  discussion  of  Figure  10, 
the  model  does  not  predict  the  experimentally  observed  decrease  in  relaxation  response  with 
increased  bolt  preload.  It  should  be  noted  that  the  measured  relaxation  is  the  total  of  all 
relaxation  mechanisms,  whereas  the  predicted  value  is  only  due  to  the  viscoelastic  relaxation  in 
the  composite  matrix  material.  Relaxation  in  the  average  tensile  normal  stresses,  az,  across  the 
pre-tension  section  in  the  bolt  was  also  observed  as  shown  in  Figure  12  but  the  corresponding 
relaxation  of  average  shear  stresses  was  negligible  for  a  20  hour  duration.  The  shift  in  the 
average  tensile  normal  stresses  is  due  to  the  application  of  static  beam  load. 

The  percentage  of  bolt  load  relaxation  (p  =  (po  -P(i))/P0  *100)  in  composite  joints  at  the  end  of  30 

hour  duration  is  shown  in  Figure  13  for  various  preloads  with  different  loading  conditions.  It  is 
observed  that  for  any  external  loading  condition,  the  bolt  load  relaxation  decreases  with 
increasing  initial  bolt  preload.  These  Findings  emphasize  the  importance  of  preload  selection.  For 
higher  preloads  (6700  and  7850  N)  the  bolt  load  relaxation  increases  with  increasing  frequency 
of  excitation,  with  lowest  relaxation  occurring  at  1  Hz  frequency.  But  for  lower  preloads  (5050 
N)  the  relaxation  decreases  with  increased  frequency  of  excitation. 

Comparing  the  percentage  bolt  load  relaxation  in  steel  and  composite  joints  for  the  duration  of 
30  hours  (Figure  14),  it  is  observed  that  only  about  1/3  of  the  total  relaxation  is  due  to 
viscoelastic  behavior  of  the  polymer  matrix  in  the  composite,  while  the  remaining  2/3  is 
apparently  due  to  the  other  relaxation  mechanisms  such  as  plasticity  and/or  slip  in  the  bolt 
threads,  which  also  occur  in  steel  joints.  Table  3  shows  some  comparisons  of  bolt  load  relaxation 
in  steel  joints  found  in  the  literature  with  the  current  predictions  for  composite  joints  extrapolated 
using  MATLAB®  [42].  It  appears  that  the  data  in  Table  3  from  references  [36]  and  [37]  are  for 
steel  lap  joints  loaded  in  dynamic  shear,  but  this  is  the  only  relaxation  data  for  dynamic  loading 
that  was  readily  available  in  the  literature.  For  the  case  of  low  cycle  loading  (i.e.  200  and  1000 
cycles),  the  current  relaxation  rates  are  far  smaller  than  those  of  in  references  [16,  17],  and  the 
current  extrapolated  relaxation  rates  only  become  signiFicant  for  millions  of  cycles  of  loading.  It 
can  be  concluded  that  the  viscoelastic  behavior  in  the  through-thickness  direction  of  composite 
bolted  joints  may  lead  to  significant  reductions  in  the  bolt  preload  over  the  service  period  of  the 
joint  depending  upon  the  design  life  and  there  is  a  need  for  long-term  experiments  on  relaxation 
in  bolted  joints. 
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CONCLUSIONS 

•  Experiments  have  been  employed  to  study  the  effects  of  various  bolt  preloads,  along  with 
the  effects  of  static  and  dynamic  external  loads  on  bolt  load  relaxation  in  composite 
bolted  joints,  and  phenomenological  models  have  been  developed  from  these  results. 

•  Finite  element  models  for  bolt  load  relaxation  in  bolted  composite  joints  based  on  a 
global/local  quasi-elastic  approach  show  reasonably  good  agreement  with  experiments 
except  that  the  experimentally  observed  decrease  in  relaxation  with  increased  bolt 
preload  is  not  predicted  by  the  models. 

•  Experiments  show  that  for  any  external  loading  condition  the  bolt  load  relaxation 
decreases  with  increasing  initial  bolt  preload,  and  these  findings  emphasize  the 
importance  of  bolt  preload  selection. 

•  If  the  bolt  preloads  are  small  enough  (as  a  percentage  of  bolt  failure  load),  applied  static 
and  dynamic  beam  loads  at  1  Hz  frequency  increase  the  magnitude  of  bolt  load 
relaxation.  However,  for  higher  bolt  preloads  the  bolt  load  relaxation  decreases  for  both 
static  and  dynamic  loads. 

•  It  is  observed  that  increasing  the  frequency  of  the  external  dynamic  load  from  1  Hz  to  5 
Hz  increases  the  rate  of  relaxation,  and  that  the  friction-induced  heating  is  at  least 
partially  responsible  for  this  effect. 

•  It  is  observed  from  the  FEA  model  that  the  average  normal  stress  across  the  pre-tension 
section  in  the  bolt  relaxes  with  time  for  the  duration  20  hours,  whereas  the  relaxation  in 
shear  stresses  is  negligible. 

•  The  FEA  predictions  of  bolt  load  relaxation  agree  well  with  the  experimental 
observations  when  subjected  to  external  static  and  dynamic  loads,  however,  more 
detailed  modeling  of  the  polymer  matrix  behavior  at  the  micromechanical  level  and 
possible  time-dependent  boundary  conditions  at  the  bolt-composite  interface  are  needed 
to  understand  the  experimentally  observed  relationship  between  bolt  preload  and  bolt 
force  relaxation. 

•  Results  of  relaxation  experiments  with  bolted  steel  joints  and  with  bolted  composite 
joints  with  Teflon  tape  thread  seal  on  the  bolt  threads  strongly  suggest  that  slip  and/or 
other  relaxation  mechanisms  in  the  bolt  threads  may  be  as  important  as  or  even  more 
important  than  the  through-the-thickness  viscoelastic  relaxation  in  the  composite  material 
being  fastened. 
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Table  1.  Constants  K  and  N  in  Equation  (1)  evaluated  using  least-squares  regression  analysis 


(a)  Composite  joint  -  Preload  only  with  no 
beam  load  (see  Figure  3) 


Preload 

K 

N 

R2 

4200  N 

0.01592 

0.3405 

0.9886 

5050  N 

0.02038 

0.1604 

0.9442 

6700  N 

0.01433 

0.2472 

0.9847 

7850  N 

0.01698 

0.1605 

0.9893 

(0  Composite  joint  -  Preload  +  250  N  dynamic 
beam  load  @5  Hz  (see  Figure  7) 


Preload 

K 

N 

R 1 

5050  N 

0.0184 

0.1293 

0.9363 

6700  N 

0.01339 

0.3027 

0.9989 

7850  N 

0.01475 

0.2271 

0.9986 

(b)  Steel  joint  -  Preload  only  with  no  beam 
load  (see  Figure  3) 


Preload 

K 

N 

R2 

5050  N 

0.01167 

0.2317 

0.9885 

6700  N 

0.007192 

0.3371 

0.996 

7850  N 

0.008088 

0.2813 

0.998 

(c)  Composite  joint  -  Preload  only  with  thread 
seal  (see  Figure  4) 


Preload 

K 

N 

R2 

5050  N 

0.03324 

0.182 

0.9667 

6700  N 

0.03397 

0.09815 

0.9019 

7850  N 

0.02452 

0.1698 

0.9771 

(d)  Composite  joint  -  Preload  +  250  N  static 
beam  load  (see  Figure  5) 


Preload 

K 

N 

R2 

4200  N 

0.01338 

0.4062 

0.9952 

5050  N 

0.02313 

0.1603 

0.9714 

6700  N 

0.01594 

0.1621 

0.9595 

7850  N 

0.009239 

0.2557 

0.9918 

(e)  Composite  joint  -  Preload  +  250  N  dynamic 
beam  load  @  1  Hz  (see  Figure  6) 


Preload 

K 

N 

R2 

4200  N 

0.02151 

0.1858 

0.9891 

5050  N 

0.02101 

0.1024 

0.9892 

6700  N 

0.01086 

0.218 

0.9564 

7850  N 

0.009212 

0.2415 

0.9345 
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Table  2.  Micromechanics  estimate  of  time-dependent  material  properties  for  unidirectional 

carbon/epoxy* 


Time 

(Hrs) 

Em(t) 

(GPa) 

E,(t) 

(GPa) 

E2(t)  =E}(t) 

(GPa) 

Gu(0  =  Gi3(t) 

(GPa) 

t>!2  (0  = 

v  U(i/0  n(0 

G23(t) 

(GPa) 

0 

3.662 

134.5 

12.07 

3.612 

0.30004 

2.535 

1 

3.607 

134.4 

11.66 

3.469 

0.30711 

2.446 

2 

3.515 

134.4 

11.55 

3.449 

0.30806 

2.434 

4 

3.495 

134.4 

11.51 

3.428 

0.30914 

2.421 

6 

3.482 

134.4 

11.48 

3.414 

0.30983 

2.412 

8 

3.473 

134.4 

11.45 

3.403 

0.31035 

2.405 

10 

3.465 

134.4 

11.40 

3.395 

0.31078 

2.400 

15 

3.450 

134.4 

11.37 

3.379 

0.31159 

2.390 

20 

3.439 

134.4 

11.34 

3.366 

0.31220 

2.382 

25 

3.430 

134.4 

11.32 

3.356 

0.31269 

2.376 

30 

3.422 

134.4 

11.55 

3.348 

0.31311 

2.371 

*  Toray  Composites  America  T800  carbon  fibers  with  assumed  viscoelastic  compliance  for  Shell  58- 
68  epoxy  matrix  material  [44] 


Table  3.  Comparison  of  rate  of  bolt  load  relaxation  due  to  dynamic  loads  from  previous  and  current 

predictions 


Number  of 
Cycles 

Plastic  deformation  of 
Threads 

(Expt.  on  steel  joints) 
(Jiang,  et  al,  [37]) 

Current 
Viscoelasticity 
Measurements  (composite 
joints)  preload  4200  N 

General  Expt. 
observation  (Steel 
joints)  (Bickford, 
f361) 

200 

41  % 

0.92% 

30% 

1000 

NA 

1.5% 

70% 

108000 

NA 

3.98  % 

NA 

3690000 

NA 

10%a 

NA 

36000000 

NA 

50%a 

NA 

a  -  Extrapolated  by  exponential  curve  fit  using  MATLAB®  [42]-  current  model 
NA  -  no  data  available 
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Figure  1.  Specimen  Configuration 
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Figure  2.  (a)  Sinusoidal  load  of  250  N  amplitude  at  frequency  of  1  Hz  (b)  Enduratec  servo-pneumatic 

testing  machine  with  3  point  bend  set-up 
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t  (hours) 

Figure  3.  Measured  bolt  load  relaxation  in  composite  and  steel  bolted  joints  normalized  to  the  initial 
preload  for  several  preloads,  along  with  predicted  (FEA)  bolt  load  relaxation  in  composite  joint  for 
preload  of  7850  N  (Note:  S  =  Steel  and  C  =  Composite) 
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Figure  4.  Measured  bolt  load  relaxation  in  composite  with  and  without  thread  seal 
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Figure  5.  Bolt  load  relaxation  for  composite  bolted  joint  with  various  preloads  and  static  250  N  beam 

load  (Fo=250  N), 
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1.0050  -I 

1  0Q00  ^  1 2.5%  (4200  N)  -0-  1 5%  (5050  N) 


Figure  6.  Bolt  load  relaxation  for  composite  bolted  joint  with  various  preloads  and  a  250  N  amplitude 
dynamic  beam  load  at  1  Hz  frequency  (F(t)=250  sin2it(l)t  N)) 
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1.005  -1 

1  - 


Figure  7.  Bolt  load  relaxation  for  composite  bolted  joint  with  various  preloads  and  a  250  N  amplitude 
dynamic  beam  load  at  5  Hz  frequency  (F(t)=250  sin2n(5)t  ,V) 
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Figure  8.  Bolt  FEA  model  (Local) 
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4220  i 


4200  ^  Preload  Only 

4180  J  — ■ — Preload+Static  125N 


Figure  9.  Predicted  bolt  load  relaxation  for  composite  bolted  joints  with  preload  4200  N 
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Figure  10.  Predicted  bolt  load  relaxation  for  composite  bolted  joints  with  various  preloads 


Po(  N) 

Figure  1 1.  Volume  averaged  macromechanical  (or  "effective")  von  Mises  stresses  in  the  immediate 

vicinity  of  the  bolt  in  composite  material 
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Figure  12.  Average  normal  stress  relaxation  in  the  pretension  section  of  bolt  when  composite  joint  is 

preloaded  to  4200  N 


12.5%  (4200  N)  15%  (5050  N)  20%  (6700  N)  23.5%  (7850  N) 


Bolt  preloads  ( Po )  as  a  percentage  of  bolt  failure  load 


Figure  13.  Measured  bolt  load  relaxation  in  composite  bolted  joint  at  30  hour  duration 
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Bolt  Preload  (P0)  as  a  percentage  of  bolt  failure  load 

Figure  14.  Bolt  load  relaxation  comparison  in  steel  and  composite  bolt  joint  for  various  preloads  under 

preload  only  condition 


53 


Part  2:  Analytical  Investigation 

INFLUENCE  OF  BOUNDARY  CONDITIONS  RELAXATION  ON  PANEL  FLUTTER 
WITH  COMPRESSIVE  IN-PLANE  LOADS 


Summary  of  Results 

The  influence  of  boundary  conditions  relaxation  on  two-dimensional  panel  flutter  is  studied  in 
the  presence  of  in-plane  loading.  The  boundary  value  problem  of  the  panel  involves  time- 
dependent  boundary  conditions  that  are  converted  into  autonomous  form  using  a  special 
coordinate  transformation.  Galerkin’s  method  is  used  to  discretize  the  panel  partial  differential 
equation  of  motion  into  six  nonlinear  ordinary  differential  equations.  The  influence  of  boundary 
conditions  relaxation  on  the  panel  modal  frequencies  and  LCO  amplitudes  in  the  time  and 
frequency  domains  is  examined  using  the  windowed  short  time  Fourier  transform  and  wavelet 
transform.  The  results  of  this  task  revealed: 

•  The  relaxation  and  system  nonlinearity  are  found  to  have  opposite  effects  on  the  time 
evolution  of  the  panel  frequency. 

•  Depending  on  the  system  damping  and  dynamic  pressure,  the  panel  frequency  can 
increase  or  decrease  with  time  as  the  boundary  conditions  approach  the  state  of  simple 
supports. 

•  The  largest  Lypunov  exponent  is  also  determined.  They  reveal  complex  dynamic 
characteristics  of  the  panel,  including  regions  of  periodic,  quasi-periodic,  and  chaotic 
motions. 

State-of-the-Art 

It  has  been  observed  that  apparently  identical  aircraft  can  exhibit  different  dynamic 
characteristics  under  the  same  flight  conditions.  This  difference  owes  its  origin  to  the  stochastic 
nature  of  structural  properties  and  the  environment.  That  is,  the  sensitivity  of  the  dynamic  system 
behavior  is  directly  linked  to  variations  in  its  physical  properties.  The  physical  properties  of 
aeroelastic  structures  are  affected  by  boundary  conditions  relaxation  and  joint  uncertainties. 
Generally,  the  main  sources  of  uncertainties  of  aerospace  structures  include: 

(v)  Randomness  in  material  properties  due  to  variations  in  material  composition. 

(vi) Randomness  in  structural  dimensions  due  to  manufacturing  variations  and  thermal 
effects. 

(vii) Randomness  in  boundary  conditions  due  to  preload  and  relaxation  variations  in 
mechanical  joints. 

(viii)  Randomness  of  external  excitations. 

The  present  work  deals  with  the  third  source  and  its  mechanisms.  There  are  many  factors  that 
affect  mechanical  joints  and  fasteners,  such  as  friction,  hardness,  finish,  and  dimensions  of  all 
parts,  and  gasket  creep  (Bickford,  1990).  Each  factor  will  vary  from  fastener  to  fastener  and  joint 
to  joint  because  of  manufacturing  or  usage  tolerances.  Moreover,  a  fastener  subjected  to 
vibration  will  not  lose  its  pre-load  immediately.  First  there  is  a  slow  loss  of  pre-load  caused  by 
various  relaxation  mechanisms.  Vibration  increases  relaxation  through  consequent  wear  and 
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hammering.  After  sufficient  pre-load  is  lost,  friction  forces  drop  below  a  critical  level  and,  if  the 
joint  is  bolted,  the  nut  actually  starts  to  back  off  and  shake  loose.  As  relaxation  occurs,  the  joint 
fails  to  mimic  ideal  boundary  conditions;  instead,  the  joint’s  properties  become  time  dependent 
and  uncertain. 

The  present  work  is  motivated  by  some  recent  results  on  the  sensitivity  and  variability  of  the 
response  of  structural  stochasticity  (see,  for  example,  Ibrahim,  1987,  and  Manohar  and  Ibrahim, 
1999)  and  by  the  recent  assessment  of  joint  uncertainties  by  Ibrahim  and  Pettit  (2004).  These 
problems  are  complex  in  nature  because  every  joint  involves  different  sources  of  uncertainty  and 
non-smooth  nonlinear  characteristics.  For  example,  the  contact  forces  are  not  ideally  plane 
because  of  manufacturing  tolerances.  Furthermore,  the  initial  forces  will  be  redistributed  non- 
uniformly  in  the  presence  of  lateral  loads.  This  is  in  addition  to  the  prying  load,  which  induces 
nonlinear  tension  in  the  bolt  and  nonlinear  compression  in  the  joint.  The  main  problems 
encountered  in  the  design  analysis  of  bolted  joints  with  parameter  uncertainties  include  random 
eigenvalues,  response  statistics,  and  probability  of  failure. 

The  combined  effect  of  uncertainty  in  the  boundary  conditions  and  spatially  variable  material 
properties  on  the  nonlinear  panel  aeroelastic  response  was  studied  by  Lindsley,  et  al.  (2002a,  b). 
It  was  shown  that  the  flutter  problem  of  aeroelastic  structures  could  be  handled  when  random 
uncertainties  are  introduced  in  the  structural  model.  Pinned  and  fixed  boundary  conditions  were 
modeled  as  limiting  cases  of  rotational  springs  on  the  boundary,  which  possessed  zero  and 
infinite  stiffness,  respectively.  Accordingly,  rotational  spring  stiffness  was  used  to  parameterize 
the  boundary  conditions.  Parametric  uncertainty  was  examined  by  modeling  variability  in 
Young’s  modulus  and  the  boundary  condition  parameter.  The  variability  in  the  boundary 
conditions  was  restricted  to  a  single  value  along  the  plate  boundary  edges  for  each  realization. 
For  values  of  the  dynamic  pressure  in  the  deterministic  limit  cycle  oscillation  (LCO)  range,  the 
variability  in  the  boundary  conditions  affected  the  plate  deflection  in  an  essentially  linear 
manner.  However,  for  values  of  dynamic  pressure  in  the  neighborhood  of  the  bifurcation  point, 
the  relationship  was  nonlinear.  Variation  in  boundary  conditions  resulted  in  a  softening  effect  of 
the  clamped  panel,  and  thus  induced  an  increase  in  the  amplitude  of  plate  oscillations. 

Structural  and  material  uncertainties  were  also  considered  in  studying  the  flutter  of  panels  and 
shells  by  Liaw  and  Yang  (1991a, b),  and  Kuttenkeuler  and  Ringertz  (1998).  For  example,  Liaw 
and  Yang  (1 991  a,b)  quantified  the  effect  of  parameter  uncertainties  on  the  reduction  of  the 
structural  reliability  and  stability  boundaries  of  initially  compressed  laminated  plates  and  shells. 
For  buckling  analysis,  the  uncertainties  include  modulus  of  elasticity,  thickness,  and  fiber 
orientation  of  individual  lamina,  as  well  as  geometric  imperfections.  For  flutter  analysis,  further 
uncertainties  such  as  mass  density,  air  density,  and  in-plane  load  were  also  considered. 
Kuttenkeuler  and  Ringertz  (1998)  performed  an  optimization  study  of  the  onset  of  flutter,  with 
respect  to  material  and  structural  uncertainties. 

A  ground  vibration  test  was  used  by  Potter  and  Lind  (2001)  to  obtain  uncertainty  models,  such  as 
natural  frequencies  and  their  associated  variations,  which  can  update  analytical  models  for  the 
purpose  of  predicting  robust  flutter  speeds.  Different  norm  approaches  were  used  to  formulate 
uncertainty  models  that  cover  the  entire  range  of  observed  variations.  It  was  found  that  the 
oo -norm  produces  the  smallest  uncertainty  and  the  least  conservative  robust  flutter  speed.  Lind 
and  Brenner  (2000)  introduced  a  tool  referred  to  as  the  “flutterometer”  for  predicting  the  onset  of 
flutter  during  a  flight  test.  The  flutterometer  computes  a  flutter  for  an  analytical  model  with 
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respect  to  an  uncertainty  description.  Brenner  (2002a)  considered  a  technique  that  identifies 
model  parameters  and  their  associated  variances  from  flight  data.  Later  Prazenica,  et  al.  (2003) 
introduced  a  technique  for  estimating  uncertainty  descriptions  based  on  a  wavelet  approach,  but 
relies  on  the  Volterra  kernels. 

The  studies  of  panel  flutter  were  concentrated  on  parametric  analysis  of  stability  boundaries  and 
the  amplitude  of  LCO  under  different  boundary  conditions.  At  the  same  time,  it  was  shown  that  a 
panel  subjected  to  a  combination  of  airflow  and  in-plane  loading  experiences  a  complex  range  of 
motions,  including  static  buckling  (divergence),  quasi-periodic  motion,  and  chaos  in  addition  to 
LCO.  Dowell  (1982)  showed  that  a  panel  under  the  combined  effect  of  fluid  flow  and  in-plane 
compression  exhibits  chaotic  motion  for  certain  values  of  some  control  parameters.  Dowell 
(1984)  observed  chaos,  via  period  doubling  and  intermittency  while  increasing  the  compressive 
in-plane  loading.  The  existence  of  multiple  attractors  and  the  coexistence  of  both  symmetric  and 
asymmetric  LCO  were  observed  by  Bolotin,  et  al.  (1998)  using  a  two  degree-of-freedom 
approximation  of  an  elastic  panel.  They  studied  the  transition  between  different  stability  regions. 
The  stability  regions  of  a  simply  supported  two-dimensional  panel  subjected  to  compressive 
loading  were  revisited  recently  by  Epureanu,  et  al.  (2004).  They  used  bifurcation  diagrams  for 
two  control  parameters  to  determine  stability  boundaries  and  Lyapunov  exponents.  The  effect  of 
damping  on  stability  boundaries  as  well  as  on  LCO  was  considered  by  Kuo,  et  al.  (1972), 
Bismark-Nasr  and  Bones  (2000),  Bolotin,  et  al.  (2002),  Pourtakdoust  and  Fazelzadeh  (2003). 
Kuo,  et  al.  (1972)  showed  that  the  edge  compression  and  viscous  structural  damping  would  result 
in  an  increase  of  flutter  amplitudes  while  the  aerodynamic  damping  would  cause  a  reduction  in 
the  flutter  amplitude. 

Relaxation  effects  in  structural  joints  cause  time-dependent  boundary  conditions  and  depend  on 
the  level  of  structural  vibration.  In  other  words,  there  are  uncertainties  in  the  boundary  conditions 
in  addition  to  a  random  field  due  to  system  parameter  uncertainties.  In  this  case,  aeroelastic 
structures  will  experience  non-stationary  time-frequency  flutter,  which  is  analyzed  using  time- 
frequency  transforms  such  as  spectrographs  and  wavelet  transform.  The  time-frequency  analysis 
techniques  have  recently  been  used  to  analyze  flight  flutter  data  by  Brenner,  et  al.  (1997), 
Johnson,  et  al.  (2002),  Staszewski  and  Cooper  (2002),  and  Yu,  et  al.  (2004).  Brenner  (1997)  used 
time-frequency  signal  representations  to  analyze  aeroelastic  flight  data.  Mastroddi  and  Bettoli 
(1999)  conducted  wavelet  analysis  in  the  neighborhood  of  a  Hopf  bifurcation  to  capture  the 
features  of  transient  responses.  In  the  neighborhood  of  aeroelastic  flutter  during  flight  tests,  the 
time  scale  decompositions  of  continuous  wavelet  transform  was  used  to  analyze  pre-  and  post- 
critical  transient  behavior  of  nonlinear  aeroelastic  structures.  Brenner  (2002b)  applied  the 
singular-value  decomposition  to  aeroelastic  pitch-plunge  wing  section  models  to  detect 
instability  and  nonlinear  dynamics  from  the  time-frequency  map. 

The  present  work  deals  with  the  nonlinear  panel  flutter  with  relaxation  in  boundary  conditions. 
The  conventional  boundary  value  problem  of  the  panel  involves  time-dependent  boundary 
conditions,  which  are  converted  to  an  autonomous  form  using  a  special  coordinate 
transformation  inspired  by  the  work  of  Qiao,  et  al.  (2000).  The  present  analysis  extends  the 
analysis  of  Ibrahim,  et  al.  (2004)  to  include  six-mode  interaction  in  the  presence  of  boundary 
condition  relaxation.  The  dynamic  characteristics  of  the  panel  and  the  influence  of  initial 
conditions  are  predicted  using  phase  plots,  FFT  plots,  bifurcation  diagrams  of  the  first  return, 
short  time  Fourier  transform,  wavelet  transform,  and  Lyapunov  exponents. 
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II.  ANALYTICAL  MODELING 
II.l  Panel  Flutter  with  Non-ideal  Boundary  Conditions 

Consider  a  two  dimensional  panel  exposed  to  supersonic  flow  as  shown  in  Figure  1.  In  order  to 
estimate  the  work  done  by  aerodynamic  loading,  the  pressure  on  the  panel  is  represented  by  the 
linear  piston  theory,  (Ashley  and  Zartarian,  1956), 


&P  =  P-P« o 


PjJl  dw  |  1  dw 
m  dx  u  dt 

_  oo 


(1) 


where  w(x,/)  is  the  panel  deflection,  which  is  a  function  of  position,  x,  and  time,  /. 
M  =Un/cin  is  the  Mach  number,  Un  is  the  undisturbed  gas  flow  speed,  ax  =  -JypJ~p1  is  the 
speed  of  sound,  px  and  px  are  the  undisturbed  free  gas  stream  pressure  and  density, 
respectively,  p  is  the  pressure  of  the  gas  flow  at  the  panel  surface,  y  =  Cp/Cy  is  the  ratio  of 
specific  heats  at  constant  pressure,  Cp ,  and  volume,  Cv . 


The  governing  nonlinear  equation  of  motion  for  the  panel  is  developed  using  Hamilton's 
principle,  which  yields  (Ibrahim,  et  al.,  1990) 
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mp  is  the  panel  mass  per  unit  area,  a  is  the  panel  length,  E  is  Young's  modulus,  h  is  the  plate 
thickness,  D  =  Eh}  /(1 2(1  -v2))  is  the  panel  stiffness,  v  is  Poisson's  ratio,  Ap0  is  the  gas 
pressure  difference  across  the  panel,  Nx0  is  the  external  in-plane  load  per  unit  span-wise  length, 
and  c  is  a  linear  viscous  damping  coefficient.  Equation  (2)  is  subject  to  the  boundary  conditions 
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where  «,(/)  and  a2(t)  measure  the  end  slopes  and  represent  torsional  stiffness  parameters  such 
that  if  a](t)  =  a2(t)  =  oo  the  panel  is  purely  clamped-clamped.  On  the  other  hand,  the  panel  is 
simply  supported  if  «,(/)  =  a2(t)  =  0 .  In  real  situations,  «,(/)  and  a2(t)  do  not  assume  these 

limiting  cases;  instead,  they  are  very  large  for  clamped  supports  or  very  small  for  simple 
supports.  In  the  dynamic  case  the  boundary  conditions  (3a,c)  are  non-autonomous.  In  order  to 
convert  these  conditions  into  an  autonomous  form,  we  introduce  the  following  transformation  of 
the  response  coordinate, 
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where  the  dimensionless  parameter  z,(t)  =  D/aa^t),  /  =  1,2,  represents  the  ratio  of  the  bending 
rigidity  to  the  torsional  stiffness  of  the  joints.  The  functions  g,(z,,z2)  and  g2(z,,z2)  are  chosen 
to  render  the  boundary  conditions  autonomous  for  the  new  coordinate  u(x,t) .  Possible 
expressions  of  these  functions  are 
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II. 2  Relaxation  of  Boundary  Stiffness 

The  relaxation  process  is  phenomenologically  modeled  based  on  experimental  results  (Bickford, 
1990).  In  this  case,  The  torsional  stiffness  parameters  are  assumed  functions  of  the  number  of 
vibration  cycles,  n  =  «(r) , 
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where  the  overbar  denotes  a  dimensionless  parameter.  An  explicit  analytical  expression  for  the 
parameters  at(n )  can  be  obtained  from  experimental  records  (Bickford,  1990),  which  reveal  a 

slow  drop  between  an  original  and  an  asymptotic  value  of  the  joint  stiffness.  An  appropriate 
elementary  function  that  emulates  this  behavior  may  be  selected  in  the  form 
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a(n)  =  a(oo)  +  [a(0)-a(oo)] 
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1  +  tanh[&Ht.  ] 


where  the  subscript  i  has  been  dropped,  and  nc  is  a  critical  number  of  cycles,  indicating  the 
location  of  the  inflection  point  with  respect  to  the  origin,  n  =  0  .  The  parameter  k  is  associated 
with  the  slope  of  the  curve  at  the  point,  n  =  nc.  The  parameters  a( 0)  and  a{ oo)  are  obtained 

from  the  experimental  curve.  The  slope  parameter  k  can  be  found  by  taking  the  derivative  of 
equation  (9)  with  respect  to  n ,  i.e., 
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One  can  write  an  expression  for  z(r)  by  using  relations  (8)  and  (10)  in  the  form 


z(r)  =  Z0ZM 


Z0-(Z0-ZJ 


l  +  tanh(-/y(r-rc)) 
l  +  tanh(^rc) 


(11) 


where  Z0  =  z(0) ,  ZK  =  z(oo) ,  %  =  — ,  and  <m>  is  the  mean  value  of  the  response 
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frequency,  which  can  be  taken  as  the  center  frequency.  The  phenomenological  representation 
given  by  equation  (11)  can  be  used  for  any  initial  preload  and  will  cause  the  panel  to  experience 
non-stationary  behavior.  Notice  that  the  relaxation  time  interval  documented  by  Bickford  (1990) 
is  very  short  for  aerospace  structural  components.  However,  the  present  analysis  serves  to  reveal 
the  dynamic  characteristics  of  panels  under  boundary  conditions  relaxation. 

II.3  Galerkin  Formulation 

Galerkin’s  method  is  applied  to  discretize  equation  (7)  by  assuming  the  general  solution  in  the 
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form  i7(T,r)  =  ^4/  „(*)<7„(r)  and  the  corresponding  weighting  functions  u(x,t)  = 

n= 1 
N 

X^»(^)^(r)  where  N  is  the  total  number  of  the  basis  functions  for  u(x,t)  ;  qn( r)  are 

n= 1 

unknown  functions  to  be  determined  (generalized  coordinates);  qn(r)  are  arbitrary  functions  of 

time  and  vF„(x)  are  the  assumed  orthonormal  mode  shapes.  The  resulting  general  differential 
equation  is 
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(13) 


The  general  solution  is  assumed  to  be  periodic  in  space, 

N 

u(x,T)  =  Yjq„(^)s\nn7rx 

n=\ 


where  N  is  the  total  number  of  modes,  qn(r)  are  the  generalized  coordinates.  It  has  been 

established  that  accurate  solution  of  the  panel  flutter  can  be  achieved  by  using  at  least  six  modes 
(see,  e.g.,  Dowell,  1966).  The  inclusion  of  six  modes  results  in  a  complicated  analysis  where 
relaxation  is  considered.  For  this  reason  we  introduce  the  simplification,  z,  =z2  =z/2,  which 

makes  the  boundary  stiffness  values  to  be  equal.  The  resulting  set  of  six  equations  may  be 
written  in  matrix  form, 


[Af(r)]{^  +  [c(C,C,A,r)]{^}  +  [jr(r,Ar0,A)]{ff}  =  [/)(r)]{^}+  X  £  {eq,q)\ 

/=I.3,5>=1JW 

i*  j*k 

+Z  Z  Z{/w*Mp} 


(14) 
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where  [A/(r)]  is  the  time  dependent  mass  matrix  and  \c(£,£,X,  r)]  is  the  damping  matrix, 

which  depends  on  the  viscous  damping  ratio  C,  ,  mass  parameter,  Q ,  and  relaxation  parameter, 

z(r).  [tf(r,./V0,A)]  is  the  stiffness  matrix,  [Z)(r)]  is  the  coefficient  matrix  of  cubic  terms,  and 

{/*}  is  the  pressure  vector,  whose  elements  are  non-zero  only  for  odd  modes  because  it  is  the 

integral  of  a  constant,  p0 ,  and  the  sinusoidal  basis  functions.  The  structure  of  these  matrices  is 

given  in  the  Appendix.  The  complete  set  of  expressions  for  all  coefficients  of  the  matrices  and 
vectors  of  equation  (4)  is  documented  in  Beloiu  (2005). 

Equations  (14)  are  solved  numerically  in  the  time  domain  for  a  typical  relaxation  curve.  The 
resulting  solution  is  given  in  terms  of  the  transformed  response,  u  ,  or  rather  in  terms  of  its 
modal  coordinates,  q , ,  /  =  1..6  .  One  should  estimate  the  modal  response  in  terms  of  its  physical 
generalized  coordinate, 

N 

w(x,t)  =  <p(x)u(x,t )  and  iv(x,r)  =  Xtfn(r)sin«/r3t  (15) 

n= 1 

where  (p  =  [x2+2 g,(z,,z2)x +g2(z,,z2)]  and  g,  are  given  by  equation  (5).  The  relationship 
between  the  physical  coordinates  q„(r)  and  the  generalized  transformed  coordinates  q„(r)  is 


™  ^  7V 

X  (r) sin  «*x  =  [x 2  +  2g,  ( z)x  +  g2  (z)]  X  (r)  sin  nnx 

n= 1 

Integrating  the  above  equation 


(16) 


n=\ 


I 


Yjqn(r)s\nn7rx 


n= 1 


dx=  J[  (x2  +2g,(z)x  +  g2(z))£qn(r)sinnxx 
L  n= 1 


dx 


gives  the  desired  relation  between  the  coordinates. 

qn^)  =  TSz)qn^) 


(17) 


where 


Uz)  =  ~ 


(2/?-l)2^2 
2  +  (2«-l)2  n2z 


n  =  1,3,5,..., 


TM  =  - 


(n-  l)2zr2 
2  +  («-l)2^-2z 


n  =  2,4,6,... 


(18) 


Therefore,  the  solution  of  equations  (14)  must  be  divided  by  Tn(z)  in  order  to  recover  the  actual 

modal  displacements.  The  next  section  presents  the  stability  analysis  and  response  characteristics 
under  different  values  of  the  dynamic  pressure  and  relaxation  parameters. 
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III.  LINEAR  ANALYSIS 

The  stability  analysis  is  carried  out  by  estimating  the  natural  frequencies  of  the  six  modes  in  the 
absence  of  system  nonlinearities  and  by  setting  the  non-homogeneous  term  in  equations  (14)  to 
zero.  The  dependence  of  the  real  and  imaginary  components  of  the  eigenvalues  on  the  dynamic 
pressure  is  shown  in  Figures  2(a)  and  (b)  for  three  different  values  of  relaxation  parameter 

(z  =  0.001,  0.1,  and  1),  damping  parameter,  ^  =  0.0,  mass  parameter  =  0.1,  and  static  axial 
load  parameter  N0  =  0.  It  is  seen  that  the  real  parts  are  zero  up  to  a  critical  value  of  the  dynamic 

pressure,  depending  on  the  value  of  the  relaxation  parameter,  z ,  above  which  one  becomes 
negative  and  the  other  positive  indicating  the  occurrence  of  panel  instability  (flutter).  Note  that 
the  value  z  =  0.0  corresponds  to  a  clamped-clamped  panel,  while  z  =  oo  corresponds  to  simple 
supports.  As  expected,  the  linear  flutter  point  decreases  for  lower  boundary  stiffness.  The 
dependence  of  the  components  of  the  first  and  second  eigenvalues  on  the  relaxation  parameter, 
z ,  is  shown  in  Figure  3  for  three  different  values  of  dynamic  pressure,  A  =  400,  450,  and  500.  It 
is  seen  that  the  eigenvalues  possess  negative  real  parts  up  to  a  critical  value  of  relaxation 
parameter,  above  which  one  eigenvalue  has  a  positive  real  part  indicating  the  occurrence  of 
flutter. 

Figures  4  and  5  show  the  boundaries  of  panel  flutter  in  terms  of  the  critical  value  of  aerodynamic 
pressure,  Acr ,  and  the  relaxation  parameter,  z  .  These  figures  depict  the  influence  of  the  in-plane 

load,  N0 ,  and  damping  ratio,  £ ,  respectively.  As  expected,  the  compression  in-plane  loading 

results  in  a  reduction  of  the  critical  flutter  speed.  The  clamped  panel  (z  « 1 )  requires  more  in¬ 
plane  compression  load  to  reach  its  flutter  speed.  With  reference  to  Figure  5,  for  all  values  of 
relaxation  parameter,  the  damping  is  non-beneficial  as  it  increases  from  very  small  values  up  to  a 
critical  value,  above  which  it  becomes  beneficial,  depending  on  the  value  of  the  relaxation 
parameter.  Figure  6  shows  the  dependence  of  flutter  speed  on  the  damping  parameter,  .  For  a 

given  relaxation  parameter,  there  is  a  critical  damping  ratio,  above  which  the  damping 
becomes  beneficial  and  the  critical  speed  increases  with  the  damping.  For  £  <  Qcr  the  damping  is 
detrimental  and  results  in  a  reduction  the  flutter  speed.  The  value  of  £cr  is  shown  by  a  small 
circle  on  each  curve  and  is  determined  by  setting  dA/d£  =  0.  The  locus  of  these  points  is  shown 
by  the  dotted  curve  in  Figure  6. 


IV.  NONLINEAR  ANALYSIS 


IV.l  Bifurcation  Analysis 

The  complete  set  of  Equations  (14)  is  solved  numerically  using  the  MATLAB©  variable  solver 
ode  15  with  relative  error  tolerance  of  KT6  and  absolute  error  tolerance  of  10~9.  The  numerical 
solution  is  carried  out  for  a  given  damping  parameter,  C, ,  and  for  different  values  of  in-plane 

load,  N0 ,  dynamic  pressure,  A ,  and  relaxation  parameter,  z  .  In  order  to  avoid  the  influence  of 

transient  motion,  only  the  last  portion  of  the  steady  state  time  history  is  taken  to  estimate  the 
state  of  the  panel.  Depending  on  the  system  parameters  and  dynamic  pressure  the  panel  may 
exhibit  different  regimes  such  as  (I)  statically  stable,  (II)  static  buckling  (divergence),  (III)  limit 
cycle  oscillations,  and  (IV)  multi-period  oscillations  and  chaos.  Figure  7  shows  these  four 
regimes  on  the  plane  of  dynamic  pressure,  A,  versus  in-plane  pressure,  -N0/tt2,  for  three 
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different  values  of  relaxation  parameter  z  =  0.001,  0.1,  and  1,  in  addition  to  the  case  of  simply 
supported  panel.  The  two  values  of  in-plane  loads  -N0ln2  =  1  and  4  represent  the  Euler 

buckling  loads  of  simply-supported  and  clamped  panels,  respectively.  It  is  seen  that  as  the 
relaxation  parameter  increases  (panel  approaches  simply  supported  case)  the  regions  (Ill)  of 
LCO  and  (IV)  multi-period  oscillations/chaos  expand.  The  dynamics  of  the  panels  along  lines  A 
and  B  shown  in  Figure  7  will  be  examined  later. 

Figure  8  shows  the  dependence  of  LCO  amplitude  on  dynamic  pressure  for  zero  in-plane 
loading  and  different  discrete  values  of  the  relaxation  parameter  z  .  Note  that,  depending  on  the 
value  of  the  relaxation  parameter,  there  is  a  critical  value  of  dynamic  pressure  at  which  LCO 
begins  in  the  form  of  supercritical  bifurcation.  The  relaxation  results  in  moving  the  bifurcation 
point  to  lower  values  of  dynamic  pressure.  Under  compression  in-plane  loading,  N0  =  -3;r2  and 

under  low  values  of  dynamic  pressure  the  panel  experiences  static  buckling  as  shown  in  Figure  9. 
As  the  dynamic  pressure  increases  the  panel  enters  a  stable  state  until  the  dynamic  pressure 
reaches  the  critical  value,  Acr,  above  which  the  panel  exhibits  LCO.  A  three-dimensional 

diagram  demonstrating  the  time  evolution  of  LCO  amplitude  and  its  dependence  on  the  dynamic 
pressure  is  shown  in  Figure  10  for  zero  in-plane  loading  and  same  parameters  as  in  the  previous 
figures. 

Under  the  relaxation  curve  shown  in  Figure  1 1(a),  the  time  history  record  of  the  total  deflection 
at  x/a  =  0.75  is  shown  in  Figure  11(b)  for  in-plane  compression  loading,  N0=-6n1,  and 

dynamic  pressure,  A  =  200 .  Over  the  whole  time  domain,  the  panel  experiences  two  different 
regimes  of  oscillations,  namely  the  growing  amplitude  LCO,  and  chaotic  oscillations. 

Chaotic  flutter  is  usually  detected  by  estimating  the  largest  value  of  the  Lypunov  exponent. 
Lyapunov  exponent  measures  the  rate  at  which  nearby  trajectories  converge  or  diverge,  and  are 
numerically  calculated  using  the  algorithm  of  Wolf,  et  al.  (1985).  Equations  (14)  may  be  re¬ 
written  in  terms  of  a  set  of  first-order  differential  equations  in  the  form 

x  =  f(x;t)  (19) 

where  x  =  {q,q}T  is  the  state-space  vector,  where  T  denotes  transpose  and  /  describes  the 
nonlinear  behavior  of  the  system.  Let  jc*(t;  jc0)  be  the  reference  solution  of  the  system  (19) 
where  *0  is  the  vector  of  initial  conditions.  In  order  to  find  the  variation  of  trajectories  in  the 
neighborhood  of  the  reference  trajectory  jc*(t),  at  each  time  step  tk  we  introduce  the 
corresponding  linearized  equation 


j  =  F(**(tk)).y  (20) 

where  F(x*(tk))  is  the  nxn  Jacobian  matrix  of  the  function  /  evaluated  at  the  reference 
solution  x  *  (tk ) .  First,  we  integrate  the  differential  equation  (19)  with  initial  conditions  x0  to 
determine  the  reference  trajectory.  Simultaneously,  the  linearized  system  (20),  with  initial 
condition  y(0),  is  integrated  for  a  small  period  of  time  A t  to  obtain  a  set  of  vectors  ^(At), 
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0=1  v»n).  The  vectors  j;(At)  are  orthonormal ized  using  the  Gram-Schmidt  procedure.  The  next 
steps  consist  of  integrating  equations  (19)  and  (20)  for  another  successive  time  intervals  A/ 
using  Jr*(tk)  and  the  orthonormal izad  vector  ^(At)  as  new  initial  conditions.  Repeating  the 

above  procedure  p  times,  the y'th  Lyapunov  exponent  can  be  obtained  as  an  average  increment  of 
variation  vector  yy(t)  during  the  test  time  A / : 


Slnr^4 

pAt£  |j>y(0;tk)|| 


(21) 


where  |||  denotes  vector  norm.  To  obtain  reliable  values  of  the  X;,  long  time  integration  is 

necessary.  In  the  present  study,  Lyapunov  exponents  are  estimated  using  a  non-dimensional  time 
increment  Ar  =  0.0001.  The  computation  starts  after  r  =  100  and  continues  up  to  r  =  1,000. 
Figure  12  shows  some  regions  of  relaxation  parameter  over  which  the  Lyaponov  exponent  is 
positive,  implying  the  existence  of  chaotic  flutter.  Note  that  extreme  positive  values  of  Lyapunov 
exponent  are  found  for  multi-period  flutter  regimes. 

The  bifurcation  diagram  shown  in  Figure  12  is  obtained  by  plotting  the  first  return  points  of  the 
panel  amplitude  for  in-plane  load  parameter N0  =  -5. 8;r2,  damping  factor  g  =  0.0001,  static 

pressure  p0  =  0,  mass  ratio  parameter  <f  =  0.1,  and  dynamic  pressure  A  =  200.  The  relaxation 
parameter  vary  between  z  =  0.0025  and  z  =  1  with  an  increment  of  Az  =  0.0025  .  It  is  seen  that 
for  relatively  small  values  of  relaxation  parameter,  z<  0.105,  the  panel  experiences  symmetric 
LCO  with  increasing  amplitude  as  the  relaxation  parameter  increases  from  the  absolute  clamped 
case,  z  =  0 .  The  figure  may  be  classified  into  the  regimes  listed  in  Table  I. 

Table  I:  Panel  flutter  regimes 


0.001  <z<0.0925 
0.0925  <z  <0.1025 
0.1025  <z<0.1700 
0.1700  <  z  <  0.1850 
0.1850  <  z  <  0.2625 
z  =  0.2650 
0.2650<z<  0.2850 
0.2850<z<  0.3775 
z  =  0.3800 
0.3800  <  z  <  0.3925 
z  =  0.3950 
0.3950  <z  <0.4025 
z  =  0.4050 
0.4050  <  z  <  0.4825 
0.4825  <  z  <  0.4950 
0.4950  <  z  <  0.5075 
0.5075  <  z  <  0.6700 
0.6700  <  z  <  0.7275 


Period-one,  symmetric  (about  z  -  axis) 

Period-one,  asymmetric  (about  z  -  axis) 

Chaos 

Period-two,  mixture  of  symmetric  and  asymmetric 
Period-one,  symmetric 
Period-five,  symmetric 
Period-one,  symmetric 

Period-one,  two,  and  three,  symmetric  and  asymmetric 

Period-eight,  symmetric 

Period-one  and  two,  symmetric 

Period-seven,  symmetric 

Period-one,  symmetric 

Period-four,  symmetric 

Period-one,  asymmetric  and  asymmetric 

Chaos 

Period-doubling,  asymmetric 
Chaos 

Period-four  and  seven,  symmetric  and  asymmetric 
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0.7275  <  z  <  0.7325 
0.7325  <z  <0.7375 
0.7375  <z  <0.8925 
0.8925  <z<l 


Chaos 

Period-six,  asymmetric 
Chaos 

Period-three,  symmetric 


For  selected  values  of  relaxation  parameter,  the  phase  plots  are  shown  in  Figure  13.  The 
corresponding  FFT  plots  are  shown  in  Figure  14  and  reveal  n  spikes  for  the  multi-period  regimes 
and  continuous  spectra  for  chaotic  motion. 


Figure  15  shows  the  bifurcation  diagram  and  the  corresponding  Lyapunov  exponent  for  X  =  250 
and  N0  =  - 6n 2 .  The  switching  from  symmetric  to  asymmetric  LCO  is  more  visible  over  the 

region  015  <z  <0.37.  After  a  window  of  symmetric  LCO,  the  motion  becomes  chaotic 
( z  >  0.6325 )  with  the  increasing  of  the  Lyapunov  exponent  as  the  relaxation  parameter 
increases.  Switching  between  symmetric  and  asymmetric  LCO,  together  with  cascades  of  quasi- 
periodic  motion,  and  chaotic  flutter  with  windows  of  periodicity  make  the  panel  behavior  very 
complex. 


Figures  12  and  15  reveal  only  a  partial  view  of  the  route  to  chaotic  flutter  during  relaxation.  To 
have  a  global  picture  over  the  parametric  space  of  dynamic  pressure  and  relaxation  parameter  for 
given  values  of  in-plane  load  and  relaxation  parameter,  the  boundaries  of  chaotic  flutter  are 
shown  in  Figure  16  for  two  values  of  N0  and  two  values  of  4".  The  dynamic  pressure  varies 
from  X  =  100  to  X  =  300  with  a  step  size  AX  =  5  ,  and  relaxation  parameter  varies  from  z  =  0  to 
z  =  l  with  step  size  Az  =  0.05  so  that  each  map  is  represented  by  41x22=902  points.  The 
Lyapunov  exponent  is  computed  for  each  set  of  parameters.  A  positive  Lyapunov  exponent 
indicates  chaotic  flutter,  which  is  labeled  by  a  black  dot.  If  all  exponents  are  negative,  the  panel 
equilibrium  position  is  stable.  If  the  largest  Lyapunov  exponent  is  zero  the  panel  experiences 
stable  limit  cycle.  Both  negative  and  zero  Lyapunov  exponents  are  labeled  by  a  blank  space. 
Figure  16(a)  shows  chaos  boundaries  for  N0  =-6n 2  and  small  damping  ratio,  C,  =0.0001.  The 

motion  is  regular  for  z<0.1.  The  relaxation  process  increases  the  chaos  occurrence  but  not 
monotonically;  instead,  a  complex  pattern  is  observed.  Higher  in-plane  loads  enlarges  the  chaos 
boundaries  as  shown  in  Figure  16(c),  however,  the  switching  of  windows  with  regular  and 
chaotic  motions  is  still  visible.  By  maintaining  the  same  in-plane  load  and  increasing  damping 
ratio  to  £  =  0.001,  the  chaos  boundaries  are  reduced  considerably,  especially  for  lower  z  as 
shown  in  Figures  16(b,d). 


With  reference  to  the  path  line  A  at  N0=  -6 n1  shown  in  Figure  7,  we  consider  the  panel 

dynamic  behavior  in  terms  of  the  bifurcation  diagram  and  Lyapunov  exponent,  shown  in  Figure 
1 7  for  relaxation  parameter,  z  =  1 .  The  aerodynamic  pressure,  X ,  is  taken  as  the  control 
parameter.  A  similar  analysis  is  found  in  Epureanu,  et  al.  (2004)  for  a  simply  supported  panel 
without  structural  damping.  The  dynamic  pressure  varies  between  X  =  1 00  and  X  =  300  with  a 
step  increment  of  AT  =  0.25 .  The  bifurcation  diagram  begins  with  the  buckled  state  of  the  panel 
up  to  X  =  109  .  The  chaotic  motion  over  the  region  109.5  <  X  <152.5  is  followed  by  a  region  of 
multi-period  oscillations  up  to  T  =  173.  Another  window  with  chaotic  motion  is  found  over  the 
range  173<T<192  followed  by  a  wide  window  with  period  three  motion.  Increasing  the 
aerodynamic  pressure,  the  motion  is  chaotic  but  the  chaos  intensity  decreases  as  suggested  by  the 
decreasing  value  of  Lyapunov  exponent.  For  T>288,  the  panel  experiences  LCO  with 
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increasing  amplitude.  Figure  18  depicts  the  bifurcation  diagram  and  corresponding  Lyapunov 
exponent  for  the  same  path  A  but  with  relaxation  parameter  z  =  0.1 ,  which  is  much  closer  to  the 
clamped  state.  The  motion  is  much  more  sensitive  to  the  change  of  control  parameter  over  the 
dynamic  pressure  range  120.5  <  A  <  156 .  With  the  reference  to  the  corresponding  Lyapunov 
exponent,  the  panel  motion  exhibits  a  cascade  of  alternating  chaotic  and  periodic  oscillations. 
The  character  of  the  motion  is  changed  abruptly  over  a  small  increment  of  A . 

With  reference  to  the  path  line  B  at  A  =  140  of  Figure  7,  we  consider  the  panel  dynamics 
behavior  by  varying  the  in-plane  load  N0.  The  bifurcation  diagram  illustrated  in  Figure  19  is 
determined  for  relaxation  parameter,  z  =  l.  It  is  seen  that  the  panel  is  stable  for  small  in-plane 
load  up  to  N0=3.3  above  which  it  experiences  LCO  over  the  range,  3.3  <N0<  3.45.  At 

N0=3A5  the  panel  experiences  secondary  bifurcation  with  symmetric  period-3  oscillations. 

Increasing  the  in-plane  load,  the  motion  switches  between  asymmetric  and  symmetric  multi¬ 
harmonic  oscillations.  Further  increase  of  the  in-plane  loading  results  in  a  chaotic  motion  with 
three  windows  of  periodicity:  4.22  <  N0  <  4.39 ,  4.77  <  N0  <  4.96 ,  and  5.58  <  N0  <  5.87  . 

For  the  same  set  of  parameters,  the  response  may  be  different  depending  on  initial  conditions. 
This  is  true  not  only  for  the  chaotic  motion  but  also  for  the  periodic  oscillations.  Figure  20  shows 
phase  portraits  for  four  different  sets  of  parameters.  Each  phase  portrait  is  drawn  for  two 
different  initial  conditions  (1)  ?,(r  =  0)  =  0.1,  q\T  =  0)  =  0,/ =  2,..,6,  and  (2)  qr,(r  =  0)  =  l, 

<7((r  =  0)  =  0,/  =  2,..,6,  respectively.  One  can  observe  from  Figures  20(a)  and  (b)  that  multi- 

periodic  oscillations  corresponding  to  initial  condition  set  (1)  becomes  a  period-3  oscillation  in 
the  case  of  initial  condition  set  (2).  In  both  cases,  the  response  is  symmetric.  For  other  sets  of 
parameters  two  asymmetric  solutions  co-exist  as  shown  in  Figures  20(c)  and  (d).  The  panel 
response  follows  one  of  these  solutions,  depending  on  the  initial  condition.  The  co-existence  of 
symmetric  and  asymmetric  LCO  is  better  observed  in  the  bifurcation  diagram  shown  in  Figure 
21  for  N0=-6n 1 ,  and  r  =  0.1.  The  response  is  asymmetric  over  the  range  233  <A  <237.5 

depending  on  the  initial  conditions;  after  that,  the  response  becomes  symmetric  independent  of 
the  initial  conditions. 

IV.2  Time-Frequency  Analysis 

The  ultimate  decision  on  whether  the  motion  is  chaotic  or  not  is  given  by  the  presence  of  a 
positive  Lyapunov  exponent.  However,  the  Lyapunov  exponent  is  characterized  by  a  slow 
convergence,  which  requires  long  time  simulations  and  large  computation  resources.  Time 
limitation  may  be  critical,  especially  when  experimental  data  is  available  for  a  limited  time 
history  record.  Therefore,  at  least  for  preliminary  investigations,  the  information  from  the 
bifurcation  diagram,  phase  plot  and  power  spectrum  is  generally  sufficient.  The  relaxation  of  the 
boundary  conditions  results  in  time  variation  of  the  panel  natural  frequencies,  and  thus  the  flutter 
becomes  non-stationary.  The  Fourier  transform  does  not  reveal  the  time  dependency  of  the 
frequency  of  panel  oscillations.  The  present  work  will  adopt  two  techniques  usually  used  for 
non-stationary  signal  analysis.  These  are  the  windowed  Fourier  transform,  known  as  the 
spectrogram,  originally  developed  by  Gabor  (1946),  and  the  Morlet  wavelet  transform.  Note  that 
both  transforms  have  time-frequency  resolution  limitations  for  the  determination  of  the 
instantaneous  frequencies.  The  windowed  Fourier  transform  relies  on  the  selected  length  of  the 
window.  Any  special  features  occur  during  short  time-scales  smaller  than  the  length  of  the 
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window,  or  with  frequencies  smaller  than  those  contained  in  the  window  are  lost  and  cannot  be 
captured  by  the  windowed  Fourier  transform.  On  the  other  hand,  the  wavelet  transform  has  the 
advantage  in  that  it  follows  the  rapid  variations  of  the  instantaneous  frequencies  since  it  adjusts 
the  length  of  the  window  according  to  the  frequency  content  of  the  signal. 

For  the  case  of  the  short  time  Fourier  transform,  a  real  and  symmetric  window  g(t)  =  g(-/)  is 
translated  by  r  and  modulated  by  the  frequency  m , 

g,.m(0  =  e'm,g(/-r),  ||g||  =  1 ,  and  \\gr,m\\  =  l  (22) 

The  short  Fourier  transform,  known  also  as  the  short  time  Fourier  transform,  of  the  panel 
deflection  q\t)  is 

oo 

Sw(t,zt)=  |  w(l)g(t  -  t  )e~'mldt  (23) 

-oo 

In  the  present  work,  the  Kaiser  (1974)  window  function  is  used.  It  has  the  following  form 


g(t)  = 


/0(/?Vi-(//r)2) 

H<r 

W) 

(24) 

0 

otherwise 

where  I0  is  the  modified  Bessel  function  of  order  zero  and  of  first  kind,  (3  is  a  parameter  that 

governs  the  shape  of  the  window,  and  T  is  the  signal  total  time.  The  spectrogram  measures  the 
energy  density  of  the  flutter  deflection  w  in  the  time-frequency  neighborhood  of  (r,c r)  given  by 

2 


Psw(j,vt)  =  |5,u'(r,cr)|2  = 


\yv(t)g{t-T)e~(a'dt 


(25) 


Figures  22(a)  and  (b)  show  two  cases  of  the  FFT  plots  and  spectrograms  of  the  panel  total 
deflection  time  history  records  for  (a)  A  =  700 ,  and  £  =  0.0001 ,  and  (b)  A  =  700 ,  and  £  =  0.02, 
respectively.  For  low  damping,  Figure  22(a)  shows  the  panel  frequency  decreases  with  time  as 
the  panel  boundary  conditions  approach  the  case  of  simple  supports.  On  the  other  hand,  as  the 
damping  increases,  the  panel  frequency  increases  with  time.  There  are  two  factors  competing 
with  each  other,  namely,  the  structure  geometric  nonlinearity  and  the  relaxation  in  the  boundary 
conditions.  For  increased  damping,  the  structure  geometric  nonlinearity  overcomes  the  influence 
of  relaxation  and  the  frequency  increases  as  shown  in  Figure  22(b). 

Alternatively,  we  will  use  the  continuous  wavelet  transform  technique  to  present  the  time  history 
records  in  two-dimensional  function  of  time  and  frequency  to  reveal  the  wavelet  modulus  and 
phase.  A  wavelet  is  a  function  with  some  special  properties.  It  should  have  a  small  concentrated 
burst  of  finite  energy  in  the  time  domain  and  exhibit  some  oscillation  in  time.  The  wavelet 
transform  can  be  regarded  as  a  cross  correlation  between  the  wavelet  and  the  panel  time  history 
record.  To  analyze  signal  structures  of  very  different  sizes,  the  wavelet  transform  decomposes 
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signals  over  dilated  and  translated  wavelets.  A  wavelet  is  a  function  ^(/)eL2(M)  with  zero 
average.  It  is  normalized  ||y/|  =  1  and  centered  in  the  neighborhood  of  t  =  0.  A  family  of 
wavelets  is  obtained  by  scaling  ^  by  s  and  translating  it  by  r  (see,  e.g.,  Mallat,  1999): 


(26) 


The  wavelet  transform  (WT)  of  the  signal  w  e  L2(R)  at  time  r  and  scale  5  is 


W’(t,s)  =  {w,v.,)  =  j  (27) 

-  \Js  V  s 

The  Morlet  wavelet  is  suitable  for  analyzing  smooth  signals.  It  is  complex,  but  not  admissible 
because  it  does  not  possess  zero  mean.  Nevertheless,  the  mean  value  is  very  small,  so  they 
generally  work  well  in  computations  even  though  they  slightly  violate  some  theoretical 
conditions.  The  mother  wavelet  for  the  Morlet  wavelet  is  given  by  the  following  function 

y/(0  =  *',/Y^lVl'|2/2  (28) 

The  modulus  of  the  wavelet  transform  is  defined  as 


K?M  =  ^(Re[»7(r,5)])2  +(lm[^;(r,5)])2 


(29) 


and  the  phase  is 


<f>(r,s )  =  tan' 


Im 

\K(r,s) 

Re 

y;(T,s)~_ 

(30) 


The  square  of  the  modulus  |^H'(r,5)|“  represents  the  energy  density  distribution  of  the  signal 


over  the  time-scale  plane,  (r,s).  On  the  other  hand,  the  phase  measures  the  relative  position  of 
the  signal  and  its  analyzing  wavelet.  The  graphical  representation  of  the  WT  modulus  in  time- 
scale  plane  is  called  scalogram.  The  modulus  and  phase  provide  the  time  evolution  of  the 
frequency  components  of  the  analyzed  signal.  Figure  23  shows  the  time  history,  scalogram,  and 
wavelet  phase  for  fixed  parameters  A  =  132,  N0=-6n 2,  and  2  =  1.  The  scalogram  shown  in 

Figure  23b  illustrates  a  large  spectrum  of  frequencies  randomly  distributed  in  time.  According  to 
Newland  (1999a,b),  the  absolute  phase  is  not  a  useful  indicator  because  it  depends  on  wavelet 
location.  Flowever,  the  rate  of  change  of  phase  with  time  in  the  same  frequency  band  is  an 
interesting  parameter  because  it  is  constant  when  the  signal  is  harmonic  of  fixed  frequency  and 
phase.  Figure  23(c)  shows  the  projection  of  phase  on  the  frequency-time  plane  and  one  can  see 
the  evolution  of  phase  with  time  does  not  maintain  a  constant  value.  A  better  visualization  of 
time-frequency  evolution  of  the  wavelet  modulus  is  illustrated  in  the  three-dimensional  plot 
(Figure  23d).  Similarly,  Figure  24  shows  the  time  history,  wavelet  modulus  and  wavelet  phase 
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for  fixed  parameters  A  =  203.5,  N0  =-6x2 ,  and  z  =  0.1.  The  wavelet  scalogram  shows  a  band 

of  frequencies  varying  about  a  dominant  component  of  co *  =  20  dimensionless  frequency.  In 
addition,  intermittent  higher  frequency  components  randomly  distributed  in  time  are  observed. 
Although  the  time  history  shows  a  certain  degree  of  repeatability,  the  motion  is  still  chaotic. 
Compared  to  the  Lyapunov  exponent,  the  wavelet  transform  is  well  suited  to  the  analysis  of  short 
duration  or  intermittent  signal  components.  However,  the  wavelet  transform  cannot  provide  a 
quantitative  tool  to  measure  chaos. 


V.  CONCLUSIONS 

The  nonlinear  flutter  of  a  two-dimensional  panel  exposed  to  supersonic  gas  flow  involving  six¬ 
mode  interaction  is  studied  in  the  presence  of  non-ideal  boundary  conditions.  The  deterministic 
study  includes  stability  analysis  in  terms  of  dynamic  pressure,  relaxation  parameter,  damping 
ratio,  and  in-plane  loading.  For  in-plane  loading  below  the  critical  buckling  value,  the  panel 
experiences  LCO  above  a  critical  aerodynamic  pressure  governed  by  the  relaxation  parameter. 
For  compressive  in-plane  loads,  the  panel  experiences  periodic,  quasi-periodic  and  chaotic 
oscillations  depending  on  the  values  of  dynamic  pressure,  relaxation  parameter  and  damping 
ratios.  Bifurcation  diagrams  of  the  first  return  and  the  associated  largest  Lyapunov  exponent  are 
estimated  by  taking  the  dynamic  pressure  or  the  relaxation  parameter  or  the  in-plane  loading  as 
control  parameters.  The  chaos  regions  represented  by  the  positive  largest  Lyapunov  exponent 
were  found  to  be  reduced  for  small  relaxation  parameter.  The  initial  conditions  were  found  to 
affect  the  behavior  of  the  panel  flutter  in  the  periodicity  and  symmetry  of  oscillations.  The  time- 
frequency  analyses  of  the  panel  flutter  was  estimated  using  the  techniques  of  spectrogram  and 
Morlet  wavelet  transform.  The  importance  of  these  transforms  is  to  reveal  the  degree  of  non- 
stationarity  of  panel  flutter  in  terms  of  frequency  variations  and  nonlinear  behavior.  The  results 
presented  herein  show  that  variations  in  boundary  fixity  clearly  affect  aeroelastic  flutter 
characteristics.  It  is  recommended  that  phenomenological  models  of  relaxation  be  integrated  into 
aeroelastic  models  of  more  complex  structures  to  ascertain  whether  relaxation  could  be  an 
important  factor  in  the  observed  variability  of  aeroelastic  behavior  between  nominally  identical 
aircraft. 
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Figure  1 .  Schematic  diagram  of  a  two-dimensional  panel  with  boundary  conditions  relaxation 
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(a) 


(b) 

Figure  2.  Dependence  of  real  and  imaginary  parts  of  the  panel  natural  frequency  on  dynamic 
pressure  for  £  =  0;^  =  0.1;  A^0  =0  (a)  real  parts,  (b)  imaginary  parts. 

— z=0.001;  z=0.1; . z=l; 
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Figure  3.  Dependence  of  real  and  imaginary  parts  of  the  first  and  second  natural  frequencies  on 
relaxation  parameter  z  for ^  =  0;£  =  0.1;No  =0  (a)  real  parts,  (b)  imaginary  parts. 
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Figure  4.  Boundaries  of  panel  flutter  on  the  plane  for  different  values  of  in-plane  load  and  for 

£  =  0.1,  ^  =  0.0001. 


Figure  5.  Boundaries  of  panel  flutter  for  different  values  of  damping  factor  showing  the  reversal 

effect  of  damping  for  £  =  0. 1 ,  jV0  =  0 . 
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Figure  6.  Boundaries  of  panel  flutter  on  the  A  plane  for  different  values  of  relaxation 

parameter,  z  ,  and  for  N0  =  0.0;  <^  =  0.1.  Dashed  curve  indicates  the  critical  damping  ratio  that 
separate  between  stabilizing  and  destabilizing  damping  effects. 


parameter  C,  =  0.0001  and  different  values  of  relaxation  parameter. 
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(I)  statically  stable,  (II)  static  buckling  (divergence),  (III)  LCO,  and  (IV)  multi-period 

oscillations  and  chaos. 
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A 

Figure  8.  Bifurcation  diagram  for  different  values  of  relaxation  parameter  for  £  =  0.0001,  £  =  0.1 , 

p0  =0,and  N0=  0. 
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Figure  9.  Bifurcation  diagram  for  different  values  of  relaxation  parameter  for  4"  =  0.0001,  £  =  0.1 , 

A,  =0,  and  N0=-3tt 2. 
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Figure  10.  Three-dimensional  plots  of  amplitudes  time  evolutions  and  their  dependence  on 
dynamic  pressure  for  £  =  0.0001,  C,  =0.1, p0  =0.0,  and  N0=0. 
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Figure  1 1.  (a)  Relaxation  of  boundary  conditions  and  (b)  time  history  record  of  panel  deflection 
at  x/a  =  0.75,  for  £  =  0.0001,  p0=  0,  £  =  0.1,  AT0  =-6tt2  ,  and  A  =  200 
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Figure  12.  Bifurcation  diagram  and  the  corresponding  largest  Lyapunov  exponent 
for  N0  =-5.$?r2 ,  C,  =0.0001,  p0=  0,  £  =  0.1,  and  X  =  200 
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Figure  13.  Phase  plots  for  £  =  0.0001,  p0=  0,  £  =  0.1,  N0  =  -5.8;r2  ,and  X  =  200  corresponding 

to  Figure  1 2  for  (a)  z  =  0.05  ;  (b)  2  =  0. 1 375  ;  (c)  z  =  0.265 ;  (d)  2  =  0.38  ;  (e)  z  =  0.4075  ; 
(f)  z  =  0.5025;  (g)  2  =  0.71;  (h)  z  =  0.8;(e)  z  =  l. 
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Figure  14.  FFT  plots  for  4"  =  0.0001,  p0=  0,  £  =  0.1,  N0  =  -5.8tt2  ,and  A  =  200  corresponding  to 

Figure  12  for  (a)  z  =  0.05;(b)  z  =  0.1375;(c)  z  =  0.265;  (d)  z  =  0.38;(e)  z  =  0.4075 ;  (0 
z  =  0.5025 ;  (g)  z  =  0.71;(h)  z  =  0.8;(e)  z  =  l. 
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Figure  15.  Bifurcation  diagram  of  the  First  return  and  largest  Lyapunov  exponent  for  £  =  0.0001 , 

po=0,  ^  =  0.1,  N0  =  -6/r2  ,and  A  =  250  . 
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Figure  1 6.  Chaos  boundaries  for  z  =  1 ,  ^  =  0.1,  p0=  0 

(a)  N0  =  -6n2,  4"  =  0.0001 

(b)  N0  =  -6n2,  4"  =  0.001 

(c)  N0  =-6.8 n2,  4"  =  0.0001 

(d)  N0  =-6.8 7T2,  4"  =  0.001 
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Figure  17.  Bifurcation  diagram  of  the  first  return  and  largest  Lyapunov  exponent  along  the  path 
line  A  of  Figure  7,  for  £  =  0.0001,  p0  =0,  ^  =  0.1,  N0  =-6/r\and  z  =  1 
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Figure  18.  Bifurcation  diagram  of  the  First  return  and  largest  Lyapunov  exponent  along  the  path 
line  A  of  Figure  7  for  ^  =  0.0001,  p0=  0,  £  =  0.\,  JV0=-6;r\and  z  =  0.1 
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Figure  19.  Bifurcation  diagram  of  the  first  return  and  Lyapunov  exponent  along  the  path  line  B 
of  Figure  7  for  4"  =  0.0001,  p0=  0,  £  =  0.1,  A  =  140, and  z  =  1 
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(c) 


(d) 


Figure  20.  Phase  diagrams  for  4"  =  0.0001,  p0=  0,  £  =  0.1,  for  two  sets  of  initial  conditions: 


(1)  q, (r  =  0)  =  0.1,  <?,(r  =  0)  =  0,/  =  2,..,6.  (2)  q](r  =  0)  =  1,  tf,(r  =  0)  =  0,/  =  2,..,6. 

(a)  Z  =  200,  N0  =  -5.8  tt2  ,  z  =  0.265 ;  (b)  A  =  200 ,  N0  =  -5.8/r2 ,  z  =  0.38 ;  (c)  A  =  \19,  N0=  -6  n1 , 

z  =  0.5 ;  (d)  Z  =  270 ,  N0  =  -6x2 ,  z  =  0.5 . 
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Figure  21.  Section  of  bifurcation  diagram  of  the  first  return  for  £  =  0.0001 ,  p0=  0,  £  =  0.1, 
iV0  = -67t2  ,  and  z  =  0.1  for  initial  conditions  (1)  qx[r  =  0)  =  0.1 ,  qt(x  =  0)  =  0,/  =  2,. .,6,  and 
initial  conditions  (2)  #,(r  =  0)  =  l,  qt(T  =  0)  =  0,/  =  2,.., 6 
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(a)  A  =  700, 4"  =  0.0001 
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(b)  A  =  700,  £  =  0.02 . 


Figure  22.  FFT  plots  and  spectrograms  for  p0  =  0,  N0  =  0,  <^  =  0.1 
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(d) 


Figure  23.  (a)  Time  history,  (b)  modulus  of  WT,  (c)  phase  of  WT,  (d)  three  dimensional  plot  of 
modulus  ofWT  for  ^  =  0.0001,  p0=  0,  £  =  0.1,  A  =  132  ,  = -6;r\  and  z  =  1 
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(b) 

Figure  24.  (a)  Time  history,  (b)  modulus  of  WT,  (c)  phase  of  WT,  (d)  three  dimensional  plot  of 
modulus  ofWT  for  £  =  0.0001,  p0=  0,  £  =  0.1,  A  =  203.5,  Na  =  -6/r2 ,  and  z  =  0.1 
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CHAPTER  3 


FLUTTER  SUPPRESSION  OF  A  PLATE-LIKE  WING 
VIA  PARAMETRIC  EXCITATION 


Summary  of  Results 

The  possibility  of  suppressing  wing  flutter  via  parametric  excitation  along  the  plane  of  highest 
rigidity  in  the  neighborhood  of  combination  resonance  is  explored.  The  nonlinear  equations  of 
motion  in  the  presence  of  incompressible  fluid  flow  are  derived  using  Hamilton’s  principle  and 
Theodorsen’s  theory  for  modeling  aerodynamic  forces.  In  the  presence  of  air  flow,  the  bending 
and  torsion  modes  possess  nearly  the  same  frequency.  Under  parametric  excitation  and  in  the 
absence  of  air  flow,  each  mode  oscillates  at  its  own  natural  frequency.  In  the  neighborhood  of 
combination  resonance,  the  nonlinear  response  is  determined  using  the  multiple  scales  method  at 
the  critical  flutter  speed  and  at  slightly  higher  airflow  speed.  The  domains  of  attraction  and 
bifurcation  diagrams  are  obtained  to  reveal  the  conditions  under  which  the  parametric  excitation 
can  provide  stabilizing  effect.  The  basins  of  attraction  for  different  values  of  excitation  amplitude 
reveal: 


•  Stabilizing  effect  that  takes  place  above  a  critical  excitation  level.  Below  that  level,  the 
response  experiences  limit  cycle  oscillations,  cascade  of  period  doubling,  and  chaos. 

•  For  flow  speed  slightly  higher  than  the  critical  flutter  speed,  the  response  experiences  a 
train  of  spikes,  known  as  “firing,"  a  term  that  is  borrowed  from  neuroscience,  followed 
by  “ refractory ”  or  recovery  effect,  up  to  an  excitation  level  above  which  the  wing  is 
stabilized. 

•  The  results  of  the  multiple  scales  method  are  verified  using  numerical  simulation  of  the 
original  nonlinear  differential  equations. 

State-of-the  Art 

Aeroelastic  structures  flutter  as  a  result  of  interaction  between  structural  and  aerodynamic  forces 
at  a  critical  airflow  speed,  which  is  referred  to  as  the  flutter  speed.  The  critical  speed  can  be 
determined  based  on  the  linear  theory  of  small  oscillations  [1-4].  However,  as  the  amplitude  of 
oscillation  increases  the  structural  forces  exhibit  nonlinear  characteristics  that  bring  the  response 
into  a  bounded  limit  cycle  oscillation  (LCO)  in  addition  to  other  nonlinear  phenomena. 
Structural  nonlinearities  can  be  geometric  that  can  affect  the  whole  structure  [5,  6]  or 
concentrated  [7].  Generally,  flutter  of  aeroelastic  structures  is  undesirable  and  the  main  goal  of 
aeroelasticians  to  suppress  and  alleviate  its  effect.  The  next  two  subsections  review  the  influence 
of  concentrated  and  geometric  nonlinearities  on  the  flutter  of  lifting  structures  exposed  to 
subsonic  airflow.  The  modeling  of  aerodynamic  loading  has  been  an  important  factor  in 
developing  the  equations  of  motion  and  will  be  outlined  in  the  third  subsection.  Related  to  the 
control  and  suppression  of  flutter  is  the  stabilizing  effect  imparted  via  parametric  excitation.  This 
interesting  phenomenon  is  assessed  in  Section  2. 

1.1  EFFECT  OF  CONCENTRATED  NONLINEARITY 

Concentrated  nonlinearity  acts  locally  in  control  mechanisms  or  connecting  parts  between  wing 
and  external  stores.  This  nonlinearity  results  from  backlash  in  linkage  elements  of  the  control 
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system,  dry  friction  in  control  cable  and  push  rod  ducts,  kinematic  limitation  of  the  control 
surface  deflection,  and  application  of  spring  tab  systems  provided  for  relieving  pilot  operation. 
Breitbach  [8]  determined  the  flutter  boundaries  for  three  different  configurations  distinguished 
by  different  types  of  nonlinearities  in  the  rudder  and  aileron  control  system  of  a  sailplane.  The 
hysteretic  damping  was  found  to  result  in  a  considerable  stabilizing  effect  and  increase  of  flutter 
speed.  Similar  effects  of  nonlinearities  due  to  friction  and  backlash  on  the  dynamic  behavior  of 
aircraft  were  reported  by  De  Ferrari,  et  al.,  [9].  The  effects  of  control  system  nonlinearities,  such 
as  actuator  force  or  deflection  limits,  on  the  performance  of  an  active  flutter  suppression  system 
were  examined  by  Peloubet,  et  al.,  [10]  and  others  [11,  12],  Reed,  et  al.,  [11]  showed  that  a 
nonlinear  system  which  is  stable  with  respect  to  small  disturbances  may  be  unstable  with  respect 
to  large  ones.  Another  important  feature  was  that  a  store  mounted  on  a  pylon  with  low  pitch 
stiffness  can  provide  a  substantial  increase  in  flutter  speed  and  reduce  the  dependency  of  flutter 
on  the  mass  and  inertia  of  stores  relative  to  that  of  stiff-mounted  stores.  A  detailed  review  of 
structural  and  aerodynamic  nonlinearities  with  more  emphasis  on  concentrated  structural 
nonlinearities  is  given  by  Lee,  et  al.,  [13], 

Free-play  nonlinearity  effects  have  been  the  subject  several  studies  [14-18].  For  example, 
Laurenson  and  Trn  [14]  investigated  the  flutter  of  a  missile  with  control  surfaces  having  free- 
play  nonlinearity.  At  a  particular  flight  speed,  the  amplitude  of  oscillation,  caused  by  external 
excitation,  starts  to  build  up.  Due  to  the  presence  of  free-play  nonlinearity  in  combination  with 
increasing  amplitude  of  oscillation,  the  effective  stiffness  of  the  system  increases  and  the  motion 
becomes  stable  at  some  limited  amplitude.  Kim  and  Lee  [16]  found  that  responses  involving 
LCO  and  chaotic  motion  are  highly  influenced  by  the  pitch-to-plunge  frequency  ratio  in  an 
airfoil  with  free-play  nonlinearity.  Experimental  studies  [17]  of  a  wing  model  with  free-play 
nonlinearity  in  pitch  showed  the  appearance  of  double  LCO.  Alighanbari  [18]  studied  three- 
degree-of-freedom  airfoil-aileron  dynamics  with  free-play  nonlinearity  in  the  aileron  hinge 
moment.  Bifurcation  analysis  indicated  various  LCO  solutions  for  velocities  well  below  the 
linear  flutter  boundary.  Depending  on  the  initial  conditions  and  air  speed,  quasi-periodic  and 
chaotic  oscillations  were  reported  for  the  aileron  motion.  In  a  series  of  papers,  Bae  and  Lee  [19] 
and  Bae,  et  al.,  [20-22]  considered  the  influence  of  structural  nonlinearities,  represented  by  free- 
play  and  bilinear,  on  various  types  of  LCO  and  periodic  motions.  The  subsonic  unsteady 
aerodynamic  forces  were  modeled  using  the  doublet-hybrid  method  originally  proposed  by  Ueda 
and  Dowell  [23].  Zhao  and  Hu  [24]  considered  similar  structural  nonlinearities  and  used 
unsteady  vortex  lattice  model  to  predict  the  LCO  of  an  airfoil  section. 

1.2  EFFECT  OF  GEOMETRIC  NONLINEARITY 

The  analysis  of  two-dimensional  airfoil  with  cubic  stiffness  nonlinearities  was  conducted  in 
references  [25-29],  Lee,  et  al.,  [25]  found  that  the  flutter  boundary  was  dependent  on  the  initial 
conditions  for  a  soft  spring,  while  for  a  hard  spring  flutter  was  independent  of  initial  conditions, 
and  both  linear  and  nonlinear  flutter  speeds  were  identical.  Further,  LCO  was  observed  for 
velocities  greater  than  the  flutter  speed.  A  jump  phenomenon  in  the  pitch  amplitude  was 
numerically  detected  by  Lee,  et  al.,  [26],  and  its  location  was  found  to  depend  on  the  given  initial 
conditions.  A  frequency  relation  was  derived  by  Liu,  et  al.,  [27]  who  observed  that  the  frequency 
and  amplitude  of  limit  cycle  oscillations  do  not  depend  on  the  choice  of  initial  conditions.  A 
secondary  bifurcation  after  the  primary  Hopf  bifurcation  was  detected  by  Liu  and  Dowell  [28], 
Furthermore,  starting  from  different  initial  conditions,  the  motion  may  jump  from  one  limit  cycle 
to  another  for  different  values  of  fluid  flow  speeds.  A  chaotic  region  was  found  by  Zhao  and 
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Young  [29]  for  certain  elastic  axis  positions,  and  the  chaotic  motion  appeared  only  at  air  flow 
speed  higher  than  the  linear  divergent  speed. 

Price,  et  al.,  [30]  studied  the  response  of  a  two-dimensional  airfoil  with  bilinear  and  cubic 
stiffness  nonlinearities.  LCO  with  period  one  was  obtained  at  velocities  well  below  the  flutter 
boundary.  In  some  cases,  where  the  airfoil  was  subjected  to  small  pre-loads,  the  motion  is 
chaotic  for  both  bilinear  and  cubic  nonlinearities.  This  was  confirmed  for  cubic  nonlinearity 
using  Lyapunov  exponents.  Sing  and  Brenner  [31]  observed  asymmetric  LCO  for  certain  values 
of  flow  speed  and  elastic  axis  location. 

A  singular  perturbation  technique  based  on  normal-form  method  was  used  to  analyze  the 
stability  of  limit  cycles  of  wing-flap  flutter  [32-34].  For  example,  Dessi  and  Mastroddi  [32] 
analyzed  limit-cycle  stability  reversal  via  singular  perturbation  and  wing-flap  flutter.  A  three- 
degree-of-freedom  aeroelastic  typical  section  with  a  trailing-edge  control  surface  was  modeled 
by  including  nonlinear  springs  in  torsional  stiffness  and  hinge  elastic  moment.  The  numerical 
analysis  revealed  the  presence  of  stable  and  unstable  LCOs,  along  with  stability  reversal  in  the 
neighborhood  of  Hopf  bifurcation.  Coller  and  Chamara  [33]  investigated  the  sub-critical  and 
supercritical  nature  of  the  flutter  Hopf  bifurcation  of  a  two-degree-of-freedom  system  with 
nonlinear  restoring  forces.  Under  certain  conditions,  the  instability  gave  rise  to  stable  LCOs 
while  for  other  conditions  unstable  periodic  orbits  emerged. 

Nonlinear  modal  interaction  of  aeroelastic  structures  can  result  in  the  occurrence  of  internal 
t  resonance  [35].  As  the  free  stream  speed  increases  the  nonlinear  solution  reveals  a  LCO  close  to 
the  initial  conditions.  For  the  flow  speed  at  which  the  aeroelastic  frequencies  are  in  3:1  internal 
resonance,  the  root-mean-square  amplitude  of  plunge  degree-of-freedom  was  found  to  increase. 
Near  1:1  internal  resonance,  the  response  grows  without  limit.  It  was  also  found  that  for  an 
aeroelastic  system  with  cubic  nonlinearity,  large  response  amplitudes  were  predicted  as  the 
system  frequencies  pass  through  a  3:1  internal  resonance. 

Marzocca,  et  al.,  [36]  considered  the  determination  of  the  sub-critical  aeroelastic  response  and 
flutter  instability  of  a  two-dimensional  wing  section.  The  analytical  model  includes  the  stiffness 
and  damping  nonlinearities  in  plunging  and  pitching  degrees  of  freedom.  In  addition  to  the 
aerodynamic  loads,  an  arbitrary  time-dependent  external  pressure  pulse  was  considered.  At  zero 
flow  speed  the  plunging  amplitudes  were  found  slightly  larger  than  those  at  the  small  flow 
speeds.  However,  this  trend  is  reversed  when  the  flow  speed  is  further  increased  and  in  such  case 
larger  amplitudes  are  experienced  near  the  flutter  speed.  For  the  flow  speed  greater  than  the 
flutter  speed  the  response  becomes  unbounded.  Recently,  Tang  and  Dowell  [37]  examined  the 
influence  of  geometric  structural  nonlinear  coupling  among  the  bending  deflection,  chord-wise 
bending  deflection,  and  twist  about  the  deformed  axis  on  flutter  speed  and  LCO  of  high-aspect- 
ratio  wings.  They  found  that  LCO  above  and  below  the  perturbation  flutter  boundary  generally 
occurs  over  a  limited  range  of  flow  speed,  depending  on  initial  conditions. 

1.3.  AERODYNAMIC  MODELING 

The  analysis  of  aeroelastic  flutter  usually  depends  on  the  methods  of  modeling  and  the  nature  of 
the  unsteady  external  disturbances  to  the  structure.  Linear  modeling  is  based  on  linear 
representation  of  aerodynamic,  elastic,  and  inertia  forces.  Nonlinear  modeling,  on  the  other  hand, 
considers  the  inherent  nonlinearity  of  any  of  these  forces  or  combination  of  them.  Ashley  [38] 
presented  the  state-of-the-art  up  to  1978  of  solved  and  unsolved  problems  of  aeroelastic  stability. 
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Resolved  problems  include  the  torsional  and  flexural-torsional  divergence  of  conventional  lifting 
surfaces,  and  subsonic  and  supersonic  flutter  of  isolated  lifting  surfaces.  It  is  believed  that  the 
unresolved  problems  include  deterministic  and  stochastic  flutter  of  nonlinear  aeroelastic 
structures  in  the  presence  of  time-varying  in-plane  excitations.  This  in  addition  to  the  stochastic 
flutter  of  linear  and  nonlinear  aeroelastic  structures  in  transonic  and  supersonic  flow  regimes. 
Recently,  some  attempts  have  been  made  to  study  the  random  response  of  typical  wing  section 
by  Poirel  and  Price  [39,  40].  The  present  work  considers  the  deterministic  flutter  of  nonlinear 
aeroelastic  lifting  structures  in  the  presence  of  time-varying  in-plane  excitations  as  a  means  to 
suppress  flutter. 

One  of  the  critical  issues  of  flutter  analysis  is  the  determination  of  aerodynamic  forces  acting  on 
an  oscillating  airfoil.  Unlike  other  vibration  problems,  these  forces  depend  on  the  airfoil  motion 
itself,  thus  resulting  in  an  interaction  between  elastic,  inertia,  and  aerodynamic  forces.  The 
magnitude  of  aerodynamic  forces  is  influenced  by  the  free  stream  Mach  number.  Several  theories 
have  been  developed  to  account  for  the  Mach  number  effect.  These  include  (a)  the  modified  strip 
analysis  [41],  (b)  the  modified  two-dimensional  loading  [42],  (c)  the  rectangular  wing  theory 
[43]  and  (d)  the  subsonic  kernel  function  [44],  The  modified  strip  analysis  is  based  on 
Theodorsen’s  work  [45],  and  is  restricted  to  non-stationary,  incompressible,  potential  flow  at 
subsonic  speeds. 

It  is  believed  that  the  early  analysis  of  wing  flutter  in  an  incompressible  subsonic  flow  was  based 
on  a  quasi-stationary  formulation.  Theodorsen  [45]  developed  a  non-stationary  formulation  that 
accounts  for  the  lag  effects  of  the  unsteady  aerodynamics  for  different  frequencies.  The  analysis 
was  based  on  the  assumption  of  harmonic  oscillations  of  the  wing.  Later,  Theodorsen  and 
Garrick  [46]  conducted  an  experimental  investigation  to  examine  the  validity  of  the  predicted 
flutter  speed  and  frequency.  Theodorsen’s  solution  was  divided  into  two  parts:  a  non-circulatory 
flow  component  that  contributes  to  the  lift  and  moment  with  a  virtual  mass,  and  a  circulatory 
component  that  takes  into  account  the  vortex  caused  by  the  wake.  The  circulation  function  is  a 
complex  function  whose  real  part  modifies  the  load  component  in-phase  with  the  angle  of  attack, 
while  its  imaginary  part  accounts  for  the  out-of-phase  load  component.  Wagner  and  Jones  (see, 
e.g.,  Fung  [3])  introduced  an  approximation  to  Theodorsen’s  function  in  order  to  predict  the 
wing  dynamic  behavior  in  unsteady  aerodynamic  field. 

Although  it  would  seem  straightforward  to  use  the  real  and  imaginary  parts  directly  in  flutter 
analysis,  this  procedure  was  found  to  give  poor  results  in  comparison  with  experimental  results 
as  reported  by  Yates  [41].  The  large  phase  angles  of  complex  circulation  functions  associated 
with  two-dimensional  flow  were  found  to  be  inappropriate  for  three-dimensional  wings.  Yates 
found  that  if  phase  angles  remained  moderately  small  the  calculated  flutter  speed  would  be 
relatively  insensitive  to  changes  in  the  magnitude  of  the  imaginary  part.  The  real  part  of  the 
circulation  function  is  usually  given  in  terms  of  the  first  eigenvalue  of  the  coupled  modes  and 
Mach  number  for  different  values  of  a  frequency  parameter  (non-dimensional  flow  speed).  The 
modified  strip  theory  includes  alterations  in  the  expressions  of  lift,  moment,  and  circulation 
function  in  order  to  approximate  the  effects  of  the  span  finite  length.  Small  disturbances  to  the 
flow  are  assumed,  which  imply  that  the  resultant  fluid  speed  differs  only  slightly  from  the  free 
stream  speed. 
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1.4  SCOPE  OF  THE  PRESENT  WORK 

The  purpose  of  the  present  study  is  to  explore  the  possibility  of  suppressing  flutter  of  a  cantilever 
wing  via  parametric  excitation  in  the  neighborhood  of  combination  parametric  resonance  of 
summed-type.  The  excitation  of  the  clamped  end  is  applied  in  the  plane  of  largest  rigidity  such 
that  the  bending  and  torsion  modes  are  cross-coupled  through  the  excitation.  This  type  of 
excitation  owes  its  origin  to  a  possible  longitudinal  vibration  of  aircraft  engines.  The  sources  of 
nonlinearities  are  due  to  the  small  in-plane  displacement  and  nonlinear  curvature.  The  paper  is 
organized  as  follows.  Section  2  deals  with  the  analytical  modeling  of  the  nonlinear  wing  under 
aerodynamic  loading.  The  aerodynamic  lift  and  moment  are  modeled  based  on  Theodorsen’s 
theory  and  the  equations  of  motion  are  developed  using  Hamilton’s  principle.  Since  the  reference 
axis  passes  through  the  elastic  axis  of  the  wing,  there  will  be  a  linear  coupling  in  the  equations  of 
motion  in  addition  to  the  nonlinear.  Thus,  the  equations  of  motion  differ  from  those  derived  by 
Ibrahim  and  Hijawi  [47]  due  to  the  presence  of  linear  coupling  and  aerodynamic  loading.  Section 
3  presents  a  linear  analysis  of  the  system  normal  modes  and  parametric  stability  boundaries  in 
the  absence  and  presence  of  airflow  at  flutter  speed.  In  the  absence  of  aerodynamic  forces,  the 
response  dynamic  behavior  is  determined  using  multiple  scales  method  in  the  neighborhood  of 
combination  resonance.  The  parametric  excitation  at  the  critical  flutter  speed  and  at  slightly 
higher  value  than  the  flutter  speed  is  considered  in  Section  4.  Section  5  presents  numerical 
simulation  of  the  original  coupled  bending  and  torsion  equations  of  motion  of  wing.  It  is 
interesting  to  find  out  that  the  parametric  excitation  results  in  stabilization  of  the  flutter  wing 
when  the  excitation  amplitude  exceeds  a  critical  value.  Parametric  excitation  of  aeroelastic 
structures  has  been  considered  for  a  limited  number  of  studies  [48,  49].  For  example, 
Lumbantobing  and  Haaker  [48]  studied  the  parametric  excitation  of  a  plunge  oscillator  and  a 
seesaw  oscillator.  Both  oscillators  are  described  by  a  nonlinear  Mathieu  equation  where  the 
origin  of  nonlinearity  is  the  dependence  on  the  flow  angle  of  attack.  It  was  found  that  by 
increasing  the  air  flow  speed  above  the  critical  value  the  equilibrium  position  is  re-stabilized. 
However,  if  the  parametric  excitation  is  due  to  longitudinal  turbulence,  the  frequency  of 
coalescence  occurs  at  a  lower  air  flow  speed  [39,  40],  Chin,  et  al.,  [49]  obtained  the  modulated 
equations  of  a  simply  supported  panel  in  a  supersonic  flow  to  calculate  the  equilibrium  solutions 
and  their  stability.  In  the  neighborhood  of  combination  parametric  resonance,  they  identified  the 
excitation  parameters  that  suppress  flutter  and  those  that  lead  to  complex  motions. 

2.  Analytical  Modeling  of  a  Plate-Like  Wing 

Consider  a  cantilever  wing  having  a  straight  elastic  axis,  z ,  perpendicular  to  the  fuselage  as 
shown  in  Figure  1(a).  The  wing  deformation  can  be  measured  by  the  bending  deflection,  u(z,t), 
and  torsional  angle,  a(z,t),  about  the  elastic  axis.  The  displacements  of  the  elastic  axis  along  y, 
and  z  axes  are  v(z,/)  and  w(z,t),  respectively.  The  deflection,  u(z,t),  is  considered  positive 
downward  and  a(z,t )  is  positive  when  the  leading  edge  is  up  (clockwise).  The  chord-wise 
distortion  will  be  neglected.  The  airfoil  is  exposed  to  an  incompressible  fluid  flow  of  speed  Um. 
Figure  1(b)  shows  the  projections  of  the  model  where  b  is  half  of  the  chord,  5,  is  the  distance 
between  the  elastic  axis  and  mid-chord,  52  is  the  distance  between  the  aerodynamic  center  and 
mid-chord,  and  83  is  the  distance  between  the  elastic  axis  and  inertia  axis. 

The  cantilever  wing  is  modeled  as  a  plate-like  beam  subjected  to  periodic  base  excitation,  T(t), 
in  the  plane  of  largest  rigidity  and  involves  nonlinear  coupling  between  bending  and  torsion 
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modes.  This  coupling  arises  mainly  from  the  fact  that  the  centers  of  mass  of  cross-sections 
undergo  a  small  but  important  displacement,  v,  in  the  plane  of  excitation.  This  displacement  is 
measured  in  terms  of  the  second-order  of  the  fundamental  bending,  u,  and  torsion,  a, 
displacements.  Under  an  initial  bending  displacement,  Am,  with  a  slight  twist,  a,  of  the  beam 
cross-section,  there  will  be  an  inevitable  very  small  displacement  Av«Am  associated  with 
these  displacements.  When  the  beam  is  released  from  this  state,  the  inertia  force  acting  through 
the  bending  displacement,  u(z,t),  will  give  a  torque  to  the  beam  cross-section.  Similarly, 
because  of  the  small  rotation  of  the  principal  planes  of  the  cross-section,  the  inertia  force 
contributes  a  bending  moment  about  the  local  plane  of  minimum  bending  stiffness  proportional 
to  the  twist  angle,  a .  The  bending-torsion  coupling  of  the  wing  can  be  realized  by  considering 
the  curvatures,  kx,  Ky,  and  Kz  about  x,  y,  and  z-axes,  [50], 


d2v  d2u  d2\  d2u 

Kx~~~d?Jr~d^<X' 


da  d\  d2u 
K.  = — — +  - 


dz  dz  dz 2 

Note  that  the  curvature  kx  is  very  small  and  can  be  set  to  zero  and  thus  one  can  write 


(1) 


(2) 


This  relationship  represents  the  projection  of  the  curvature  of  the  beam  in  the  plane  Oxz  on  the 
Oy-axis,  which  gives  the  curvature  in  the  plane  Oyz,  see  Figure  1.  It  is  clear  that  the 
displacement  v  can  be  expressed  in  terms  of  u  and  a  through  double  integration  of  relation  (2). 
Furthermore,  one  can  approximate  the  curvature  about  z-axis  by  the  linear  relationship 
k.  =  -da/dz .  Neglecting  the  extension  of  the  beam  elastic  axis,  the  lateral  displacement,  v(z,r), 
and  axial  drop,  w(z,/),  may  be  expressed  in  terms  of  the  bending  deflection  and  the  twist  angle 
in  the  form: 


v(z't>= 

The  kinetic  energy  is: 


T 


d\  dY 
~dt+lh 


dxdydz  +  \\I\^]  dz 


(3) 

(4) 


(5) 


where,  p  is  the  wing  density,  70  is  the  mass  moment  of  the  inertia  (about  the  inertia  axis)  per 
unit  length  and  /  is  the  wing  length.  The  strain  potential  energy  due  to  bending  and  torsion  is 
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v  = 


dz  +—cGJ 
2 


dz 


(6) 


where  E  is  Young's  modulus,  /  is  the  area  moment  of  inertia  of  the  wing  cross-section  about  y 
axis,  J  is  the  polar  moment  of  inertia  of  the  wing  cross-section  about  z  axis,  G  is  the  modulus 
of  rigidity,  and  c  is  a  correction  constant  due  to  the  noncircular  cross-section  of  the  wing.  The 
curvature  is  expressed  up  to  cubic  order. 


For  incompressible  flow,  the  aerodynamic  lift  per  unit  span  is  obtained  based  on  Theodorsen’s 
theory  [45]  in  the  form 
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The  aerodynamic  moment  about  the  elastic  axis,  z ,  per  unit  span  is 
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where  pw  is  the  air  density,  a  =  8l/b ,  B(k)  is  the  circulation  function,  which  depends  on  the 

reduced  frequency  parameter,  k  =  b(o/UK,  where  co  is  the  natural  frequency  of  the  wing  coupled 

modes  and  will  be  determined  later.  The  circulation  function  is  a  complex  quantity  represented 
by  [3,51]: 


B(k)  =  F(k)  +  iG(k) 


(9) 


where  F(k)  =  1  - 


0.165A:2 
k2  +0.00207 


given  by  Jones  [51]. 
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Applying  Hamilton's  principle 
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and  carrying  out  the  variational  process,  gives  the  equations  of  motion  and  associated  boundary 
conditions.  The  boundary  conditions  are 
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«lr=0  = 


da  _da  du 

dz  lz=0  dz  Z=l  °’ 


d2“l  =  =0 

dz2  2=1  dz3  z=l 


(11) 


Considering  only  the  first  mode  in  bending  and  torsion,  expanding  the  solutions  in  terms  of  the 
generalized  coordinates  and  mode  shapes 


u(z,t)  =  u0(t)f(z); 


a(z,0  =  a0(0<t>(z) 


(12a, b) 


where  /(z)  =  cosh(l.875z//)-cos(l.875z//)-0.734[sinh(l.875z//)-sin(l.875z//)],  and 
<})  (z)  =  sin  (txz  /  2/)  . 

Applying  Galerkin’s  method  gives  the  following  two  coupled  nonlinear  ordinary  differential 
equations 
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c7  (\  +  2a)nb2pxB(k)U„u0  (t)  +  2cAmu0  (/) d0  (/)w0  (t)  +  c5mu0  (t)Y(t)  =  0  (14) 

where  m  is  the  wing  mass  per  unit  length,  Ia=I0  +  mbl ,  Ia  is  the  wing  mass  moment  of  the 
inertia  (about  the  elastic  axis)  per  unit  length,  and  Sa=md3,  Ku  =  12.3596£/,  //4 , 
Ka=n2cGJ/4l2,  Ku3  =  80.$579EIy / 16 ,  c,  =4.597,  c2  =0.222567,  c3=0.42,  c4  =0.4552, 
c5  =0.84129,  c6  =0.677861,  and  c7  =1.3557. 

Rearranging  the  equations  of  motion,  introducing  the  non-dimensional  parameters  u=u0/b, 
Y  =  Y/b,  t  =  co at ,  and  linear  viscous  damping  with  damping  factors  t^u  and  CSa,  equations  (13) 
and  (14)  take  the  non-dimensional  form 
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where  a  prime  denotes  differentiation  with  respect  to  the  non-dimensional  time  t,  the  subscript 
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0  is  dropped,  £a  = 
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Equations  (15)  and  (16)  will  be  solved  for  linear  modal  analysis,  parametric  stability  boundaries, 
and  response  LCO  in  the  absence  and  presence  of  air  flow. 

3.  MODAL  ANALYSIS  AND  PARAMETRIC  STABILITY 


The  purpose  of  the  linear  analysis  is  to  determine  the  critical  flutter  speed  and  the  corresponding 
coupled  frequencies  on  the  system  parameters  and  airflow  speed.  Under  parametric  excitation  it 
is  important  to  determine  stability  boundaries  in  the  neighborhood  of  combination  parametric 
resonance  at  the  flutter  speed.  The  linear  equations  of  motion  are: 
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In  the  absence  of  parametric  excitation,  equations  (17)  and  (18)  are  linearly  coupled  and  may  be 
solved  for  the  eigenvalues  for  the  following  parameters:  air  mass  parameter  p  =  l/38,  elastic 

axis  location  a  =  -0.36,  inertia  axis  location  xa=  0.024,  inertia  ratio  ra=0.62,  bending-to- 
torsion  frequency  ratio  r  =  0.5,  and  damping  factors  C,a  =0.001,  C,u  =0.001.  To  estimate  the 
critical  flutter  speed,  the  solution  of  equations  (17)  and  (18)  is  expressed  in  the  form: 
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u  =  u„e'm  and  a  =  <x„e" 


(19) 


where  to  =  co/toa  is  the  nondimensional  natural  frequency  of  the  wing  coupled  modes. 

Substituting  equations  (19)  into  equations  (17)  and  (18),  and  setting  the  parametric  excitation  to 
zero,  gives 


-(l  +  p)co2+2/ 

l  K 


CD  +  A*2 


+ 


|"c6  (*a  -a\i)(£>2 +icb  -^[l  +  flU)(l  -2a)]©+  2c6  ^^|ao  =  0 


(20) 


-c7(xa  -op)co2  -/'c7  (l  +  2a)-^Z?(£)a>  ua  +  |-l--^-^  +  o2  jco2  + 


2^-|  j  "2 a 


KK 


+  \--a 


2 

a'  a 


kr 


co  + 


t  ,,  ,  \V-BUc) 


klfi 


>a„  =0 


(21) 


Substituting  for  B(k )  from  the  equation  (9),  replacing  k  =  ^SL—  =  k  to,  and  rearranging, 

(/,  to„ 

OO  d 

gives  the  characteristic  equation: 


AA(ka,  to)  BB(ka,u) 

DD(ka ,  to)  EE(ka,(o) 

where  the  elements  of  the  determinant  are  given  in  Appendix  A.  Traditionally,  (see  e.g.,  Fung 
[3]),  the  flutter  speed  and  flutter  frequencies  are  determined  by  setting  the  real  and  imaginary 
parts  of  the  determinant  (22)  to  zero.  By  giving  a  series  of  values  of  k  the  corresponding  values 
of  a>  were  determined  from  the  real  and  imaginary  parts  of  the  determinant.  The  intersection  of 
the  curves  corresponding  to  the  two  parts  establishes  the  flutter  frequency  and  flutter  speed. 
Instead,  we  will  numerically  solve  the  determinant  (22)  as  a  function  of  the  flow  speed.  In  view 
of  the  presence  of  a)  up  to  forth  order  in  the  denominators  of  circulation  functions  F  and  G , 
the  resulting  frequency  equation  is  found  to  be  of  order  12.  Careful  inspection  of  this  equation 
reveals  that  the  free  terms  (coefficient  of  co°)  and  coefficients  of  lower  order  terms  co,  co1  up  to 
fourth  or  fifth-order  are  extremely  small,  resulting  in  almost  zero  roots.  The  numerical  solution 
gives  two  different  complex  roots  and  their  conjugates.  The  values  of  the  remaining  roots  were 
found  to  be  almost  near  zero.  The  real  components  («„ /co^  and  (ro„/(oa)2  represent 
dimensionsless  natural  frequencies  of  the  coupled  modes,  and  the  imaginary  components 
(^m)i 1  and  (^/m)2  ^  “a  stanc*  f°r  the  corresponding  damping  factors.  The  dependence  of  the 
eigenvalue  components  on  the  flow  speed  parameter  Ux  /bwa  is  shown  in  Figure  2(a)  when  both 
circulation  functions  F  and  G  are  considered.  It  is  seen  that  -(^,m),/coa  is  always  negative 


=  0 


(22) 
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while  -(^m)2/coa  switches  to  a  positive  value  at  the  critical  airflow  speed  U„/bo)a  =5.04 
indicating  the  occurrence  of  flutter. 


If  one  ignores  the  imaginary  component  of  the  circulation  function,  i.e.,  G  =  0,  the  resulting 
natural  frequencies  shown  in  Figure  2(b)  are  found  very  close  to  those  shown  in  Figure  2(a).  For 
the  case  of  Figure  2(a)  the  critical  flutter  speed  ratio  is  ((/„ /be oa )  =  5.04 ,  while  for  the  case 

of  Figure  2(b)  (UK  / 6coa)G=0  =  4.965  .  The  error  between  the  two  flutter  speeds  obtained  in  both 
cases  is  1.5%.  Since  this  error  is  less  than  2%  we  will  ignore  the  contribution  of  G  . 

It  is  interesting  to  note  that  at  zero  flow  speed,  the  two  natural  frequencies  possess  the  ratio  0.5 . 


Under  parametric  excitation,  Y  =  T0cos(Qx),  where  Q  =  Q/a>a,  equations  (17)  and  (18)  are 
considered  for  stability  analysis.  Consider  the  solution  of  equations  (17)  and  (18)  in  the  form 


,  / 

'nx] 

u  =  ax  sin 

+  b.  cos 

l 2  J 

1  1 

.  2  J 

a  =  A.  sin 


(Qx)  (Qx 

—  +5,  cos  — 


2  ) 


V  2 


(23a,  b) 


Substituting  equations  (23)  into  equations  (17)  and  (18)  and  equating  the  coefficients  of 
sin(f2x/2)  and  cos(Qx/2)  gives  a  set  of  four  linear  homogeneous  algebraic  equations.  The 
determinant  of  the  coefficients  of  these  equations  establishes  the  frequency  equation: 


Q 


B(k)ix 


c7  Q2  /  x  1  ,,  - 

•— (xa-ap)  +  -c5T0Q 

_c_1mot(l+2a) 

2  k„ 


^(l  +  B(*)(l-2a)) 

O2  r„2-i^-^i(l  +  2a)-^-ia2QV 

4  k2a  32  4  h 

^(\-2a-B(k)(\-4a2))  +  r>nt;a 


Q 


B(k)ii 


+  r^ 


/ 


^Si(l  +  2  a) 

2*« 

■^-(*«  ~av)-^c5Y0n2 


^^(\  +  B(k)(\-2a)) 


2k 


cbn 


4  (*a  ~aV)~\ C3YqQ2  +  - 


-( 0-: 2« -  «<*>' 0  - )). +  , r„’l fi<;„ 


mn(l+2a)_a!jx 


r  - 


— -a2Q2\i 
32  4 


=  0 


(24) 


4 


4 
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The  roots  of  this  equation  depend  on  the  excitation  amplitude  and  flow  speed.  These  roots  are 
known  to  be  the  boundary  frequencies  of  the  instability  regions.  The  instability  regions  are 
bounded  by  periodic  solutions  of  period  2 T ,  where  T  is  the  period  of  parametric  excitation.  The 
stability  boundaries  belonging  to  zero  flow  speed,  UK/ba>a  =0,  are  given  by  the  solid  curves  in 

Figure  3,  while  those  belonging  to  the  critical  flow  speed,  Um /b(oa  =  4.965 ,  are  shown  by  the 
dashed  curves.  At  the  critical  airflow  speed,  the  stability  boundary  touches  the  frequency  axis  at 
excitation  frequency  ratio  Q/(coa  +cou)  =  1.04  regardless  of  the  damping  of  the  structure.  In  the 

absence  of  airflow,  there  are  two  instability  regions,  a  narrow  region  at  excitation  frequency  ratio 
less  than  1  and  a  wider  one  at  excitation  frequency  greater  than  one.  At  and  above  the  critical 
speed  there  is  only  one  instability  region  and  it  is  wider  than  those  below  the  critical  speed.  The 
bottom  of  each  instability  regions  at  zero  air  speed  moves  away  from  the  frequency  axis.  Inside 
instability  regions,  one  should  consider  the  system  inherent  nonlinearities  that  are  responsible  to 
bring  the  response  amplitude  into  a  bounded  LCO. 

4.  NONLINEAR  ANALYSIS 


4.1  RESONANCE  CONDITIONS 

The  nonlinear  modal  interaction  is  examined  by  considering  the  coupled  nonlinear  equations  of 
motion  (15)  and  (16)  using  multiple  scales  method  [52].  The  equations  of  motion  are  first  written 
in  terms  of  the  linearized  principal  coordinates.  Consider  the  equations  of  motion  in  non- 
dimensional  form: 


where  u(x)  = 


f«(T)l 

laCOj 


F  = 

’  Fnl 


M  (u(t))  u'(T)  +  K(*)u(t)  =  F„  (u(x),  u'(T),  Y\x),  k) 

( fnl 


(25) 


u 


n2 


-,  r(T)  =  rosin(Qx),  Q ~ , 


M  = 


2 


72—2 


c7(xa-an)  +  cAua  p 


c<,(x 

:a-op)-c2«a 

(r>  2c/S<*>  ) 

,  p 

f  1  2)  c4  -2 

,  K  = 

a 

+7 

a 

-  +  a  +-7M 

U  J  ra2  J 

0  1  (1  +  20)  **  «(*) 

^  Ka'a  / 

(  \ 

rS.+f-Bd)  u-^[l  +  BU)(l-2a)]a'-c,d!u‘u+2c2u'a’a-c,d‘r‘u‘-c,ar- 


kr 


1  » 


k  r2 

Kara 


£(£)  +  !  --a 


Kr« 


a' -l^rWa'u  -^uY" 
r  r 


Multiplying  equation  (25)  by  the  inverse  of  mass  matrix,  gives 


Iu'  +  M~'Ku  =  M"'Fnl  (26) 

where  I  is  the  identity  matrix.  The  inverse  of  the  mass  matrix  is  then  expanded  into  a  Taylor 
series  up  to  cubic  order.  The  normal  mode  natural  frequencies  are  determined  from  the 
conservative  system,  Iu"  +  Ku  =  0 ,  and  are  given  by: 
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®?(*)  = 


_ P3(*)- >//><(*) _ 

2[8ra  (1+H)  +  (1  +  8°2)(1  +  M)H-8C6C7(:,t:a 


(27a) 


(022(k)  = 


_ _ 

2  [8ra  0  +  p)  + 11  + 8tf 2  )(1  +  p) P -  Kci  (xa  - «m)2  ] 


(27b) 


where  p2(k)  and  p4(k)  are  given  in  Appendix  B.  ©,=co,/(Da,  and  co2=co2/oc>a.  The 
corresponding  modal  matrix  is: 


P(*)  = 


1 

KP\(k) 


1  N 

Pitt); 


(28) 


where  p^k)  and  p2(k)  are  the  modal  ratios  given  in  Appendix  B.  The  following  coordinate 
transformation  is  introduced: 


u(t)  =  P(*)v(x) 


(29) 


v,(x)] 


where  v(t)  =  |v*  ^  j  are  the  principal  coordinates.  The  equations  of  motion  become: 


Iv’(t)  +  A(*)v(x)  =  Fnl  ( V|  (t),v2  (x),vj  (x),v2 


(30) 


where  A (k)  = 


(02(k)  0  ^ 

0  ©2^) 


.  Fn.(vpv2,v;y2,^)=  y0u  =  y0q2. 


nil J 


=e{fev2  +^4v;2  +95v'22  +^6v!v2  +  W\v2  +%v'v2)v,  +  (q9\'2  +  ql0\[\'2  +  quV2)\2  +ql2V] 

+<7l3V2  +  ?14VI  +  9l5V2  +  (‘7l6V2  +  tfl7V'l  +  ?18V2  )  V!  +  (?19V!  +^20V2)V2  “ (^lV2  +  ?2VI  )  Y0u  SM  (fix)} 

Kl2  =e{(-S3V2  +54V!2  +^5V22  +  *6V1V2  +  J7v|v2  +  Sgv'2V2  )  V,  +  ( S9\'2  +  5,^  +  JUV?  )  V,  +  J^vj 
+J13V2  +  514VI  +  S\5V2  +  (■S16V2  +  Sn  VI  +  ■*! 8V2  )  Vl  +  (*19V!  +  S20V2  )  V2  -(51  V2  +  S2  V1  )  Y0u  S'n  (^X)} 


and  a  prime  denotes  differentiation  with  respect  to  the  non-dimensional  time  x.  qr..q20,  and 
5,...520  are  coefficients  that  depend  on  the  flow  speed.  The  solution  of  these  equations  is 
expressed  in  terms  of  two  time  scales  T0  =  x  and  Tx  =  sx 


vl(t,s)  =  v10(r0,7;)+8v11(r0,7;)+...,  v2(x,s)  =  v20(r0,7;)+8v21(7’0,7;)+...  (3ia,b) 
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Substituting  equations  (31)  into  equations  (30)  and  collecting  terms  of  the  same  order  of  8,  the 
following  equations  are  obtained: 

Equations  of  order  e°  are: 


£>nVrn+C0,2V,n  =0 


'0  v10 


10 


Do  V20  ^2^20  —  ^ 


Equations  of  order  8  are: 


A>v„  +(&i  vM  —  (^3  v20  -\-q^D0vl0  +q5  D0v20  +  q6D0\]0D0v20  +  <77£\)V,0v20  +  *7sA)v2ov2o)  vio 

+  (^9Z)0Vio  +  (j[\oD0\ iqDqV 20  +  *7l  1^0^20  )  ^20  1 2 Dq V i o  +  9l4V10  #13^0^20 

^”(^16^20  ^~9l7^\)^10  ^18^0^20  )  ^10  +  (*7l9^0^10  ^20^0^20  )  ^20 

-^'(Wio+wMu  (e~inT°  -e'nr°  )-2D0Z),v,0  j 

^0V21  V2I  —  (,S3V20  "*‘,S4^0VI0  "*'‘S5'A>V20  ,y6^0V10  A)V20  'S7-^0VI0V20  ,y8'^0V20V20  )  VI0 

+  ( 59  D0  V  |  o  +  J|0Z)qV|0Z)0V20  +  ^iiT^qVjo)  Vjq  +5|2i?0V|0  +  S|4V,0  +  +  S,jV2( 

(‘yi6V20  +  S\?DoV\o  *^1 8 2-^0 V 20  )  V 1 0  ( ^1 9 ^0 V ] 0  S 20^0^ 20  )  V 


4 


~2/(5'v2o+s2vio)Io,,  (e"nT° -e'm°)-2 D0Z),v 


10 


(32a) 

(32b) 


(33a) 


(33b) 


The  zeroth-order  solution  is: 


< 


vl0(T0Jl)  =  A[T,]e'^+A[Tl]e-r^-,  y20(T0,Tl)  =  Q[Tl]ei^  +Q[Ti]e~m'T°  (34a, b) 

Equations  (33)  contain  secular  terms,  which  give  rise  to  following  resonance  conditions: 

(a)  Combination  parametric  resonance  Q  =  co,  +  co2 ,  i.e.,  Q  =  co,  +co2 

(b)  One-to-one  internal  resonance  co,  =  (02 

(c)  One-to-three  internal  resonance  co,  =  3co2 

(d)  Three-to-one  internal  resonance  3co,  =  co2 

The  third  and  fourth  resonance  conditions  are  excluded  because  Theodorsen’s  theory  is  only 
applicable  near  flutter  speed,  i.e.,  when  co,  »  co2 .  The  first  two  cases  will  be  considered. 

4.2  RESPONSE  AT  ZERO  FLOW  SPEED 

The  results  of  the  previous  section  revealed  that  the  normal  mode  frequencies  in  the  absence  of 
air  flow  are  not  equal  and  possess  the  ratio  co,  /co2  =  0.5 .  Setting  secular  terms  corresponding  to 

the  combination  parametric  resonance,  Q  =  co,  +  to2 ,  to  zero  and  introducing  the  detuning 
parameter  cr. ,  defined  by  Q  =  co,  +S2  +8(7, ,  gives  the  solvability  conditions 
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1  .  ^ 


~  icl\QY^'  +  ty|2<M  +  (3^.4  +  iq„  CO,  +  q,  cof )  ^2/l  +  2(^3  +  /^19co,  +  qr5(o2  )AQQ-  2m, A'  =  0 


(35a) 


^  is2AY0ue'aJ'  +  /-SI3(O20  +  (3^5  +  /520(02  +  ^,(02 )  02£?  +  2  (s,6  +  w18co2  +  s9(o2  )  QAA  -  2m2Q'  =  0 

(35b) 

Introducing  the  following  polar  transformations 

A  =  V*'  i  Q  =  b2e*>  (36a, b) 

and  separating  imaginary  and  real  parts  gives  the  following  first-order  differential  equations: 

2b,  —  qi2b,  +  q,jb,  +  2q,9b2b,+-^^  b2Y0u cos((T(,7]  —  p,  —  p2)  (37a) 

263, 6,p;  =  -(3<7I4  +  ?4632  )b,  -2 [q3  +  q5(b22 ] b,b\  +  Sl b2Y0u  sin (ae7J  -  P,  -  P2 )  (37b) 

2b2  =  snb2  +  s20b2  +2snb,  b2  +  b,Y0u  cos(a,7j  — p,  — P2)  (37c) 

2cd262P2  —  —  2(^,5  +59(o,  }b2b,  —  (3s,5  +  ^|,(o2  )^2  ^  ~b,Y0u  sin  (cre7|  —  p,  —  P2 )  (37d) 

The  above  four  algebraic  equations  contain  three  unknowns  and  the  dependence  on  the  detuning 
parameter  will  not  appear  explicitly.  In  this  case,  it  is  possible  to  combine  two  of  these  equations 
into  one  through  the  following  transformation: 

ye = -Pi -p2  and  pi+p^-y:  (38) 

where  a  prime  denotes  differentiation  with  respect  to  T, .  The  resulting  equations  are: 


2b[  =  qnb,  +  qt7b2+2qt9b2b,  +^±-b2Y0u  cos(ye) 

2to, 

2b2  =  sl3b2  +  s20b2  +  2s,gZ>,  b2  +  2  b ,  Y0u  cos  (ye ) 

Z(02 


(39a) 

(39b) 


y,  =  ct  - 

l  e  e 


q,b2+s2  b, 


4(0,  b ,  4co2  b2 


You  sinye  + 


—  2 


39,4+94(0,  .  ^,6+59C0| 


2co 


1 

—2 


CO, 


35l5+5„(02  1  g3+g5  C02 
2(0-,  63. 


(39c) 


Equations  (39)  are  numerically  integrated  for  different  values  of  excitation  amplitude  in  the 
neighborhood  of  combination  parametric  resonance,  i.e.,  Q  « co,  +  co2 ,  internal  frequency  ratio 

co,  /co2  =  0.5  ,  and  the  external  excitation  detuning  parameter  <Je  =  0.  Over  a  very  small  range  of 
excitation  amplitude,  0  <  F0  <  0.0008 ,  the  zero  equilibrium  is  dynamically  stable  and  the 

damping  energy  overcomes  input  energy.  At  a  threshold  value  of  excitation  amplitude  the 
equilibrium  position  loses  its  stability  and  both  bending  and  torsion  amplitudes  achieve  steady 
state  values  over  the  excitation  amplitude  range  0.0008  <  Y0  <  0.001 8 .  Above  that  range  and  over 

the  region  0.0018  <  }q  <  0.003,  the  response  amplitudes  experience  Hopf  bifurcation.  The 
response  then  experiences  another  bifurcation  in  the  form  of  period  doubling  over  the  region, 
0.003  <  Y0  <  0.0055 .  Above  that  region,  the  response  exhibits  regular  periodic  (but  not 
harmonic)  oscillations.  Figure  4  shows  the  bifurcation  diagram  using  the  initial  conditions 
«(0)  =  0.1,  a(0)  =  0.1,  z7'(0)  =  0,  a'(0)  =  0,  and  system  parameters  p  =  l/38,  a  -  -0.36 , 

xa  =  0.024 ,  ra  =  0.62 ,  r  =  0.5 ,  C,a=  0.001 ,  C,u  =  0.001 ,  and  d  =  0.1 .  The  bifurcation  diagram 

reveals  the  upper  and  lower  pairs  of  local  maxima  and  minima  by  taking  the  excitation  amplitude 
as  the  control  parameter. 

4.3  RESPONSE  AT  CRITICAL  AND  POST-CRITICAL  FLUTTER  SPEEDS 

Introducing  the  external,  <re,  and  internal,  <r,,  detuning  parameters,  defined  such  that 

Q  =  co,  +  co2  +  e<je  and  co2  =  to,  +  scr, ,  gives  the  solvability  conditions 


+(<J,6  +i<lA +  - iq +iqsa2  +9,0w,w2  -q5a22)e2io‘T'Q2A 

+  (3<7i5  +  +  <7|  1®2 )  e<1'  Q2Q  +  2  (qt  6  +  q9l of  +  iqXi  c52 )  e,'s‘r'  QAA 

(40a) 


Introducing  the  following  polar  transformations 


Q  =  b2e ,p! 


(4 1  a,b) 


and  separating  imaginary  and  real  parts  gives  four  first-order  differential  equations.  By  setting 
the  left-hand  sides  of  these  equations  to  zero,  one  obtains  four  algebraic  equations  containing 
three  unknowns,  since  the  two  phase  angles  are  combined  as  one  unknown.  In  addition,  the 
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dependency  on  the  detuning  parameter  does  not  explicitly  appear.  In  this  case,  one  may  introduce 
the  following  transformation: 

ye=^-P,-P25  y,=cy(7;-Pi+p2,  p; +P2  =<re p;-p;=a,-y;  (42) 

where  a  prime  denotes  differentiation  with  respect  to  Tx .  The  resulting  equations  are: 


2b[  -  qx2bx  +  qx7bx  +2qx9b2bx  +^zrb2Y0u  cosyc  +  b2  cosy,  +-^L[qr2oco2  cosy,  +C,  sin  y, ] 

^COi  0)i 


il3®2 


CO, 


b}b. 


+-^[C2cosY/  +Qsiny,]  cos(ye +y,)  +C46,622cos2y,+C5^22sin2y,  (43a) 


CO 


2b'i  =  sX3b2  +s20bl  +2 sxsbxb2  +z^=rbxY0u  cos(ye)  +  ^^-bx  cosy,  +-^[.5,7(0,  cosy,  -C6  sin y,l 

ZCD,  CD,  CD, 


S^CD, 


b^b  b 

+-^-l[C7  cosy,  -C8  sin  Y,]+^3-^T0a  cos(ye  -y,)  -b2b2  (Cg  cos2y,  +CI0  sin  2y, )  (43b) 


2co, 


r'e=cre  +  A5j!dbb  {~2q"(0*b* sin  r'  +  b'  [C< cos r‘  +  2s^  sin ]  +2snAb!  sin  y, 

+bt  (cnA  cos  y,  -  2 q20al  sin  y, )  +Cl3bxb2  cos  y,  +CX4bfb2  sin  y,  +  2 CX5tfb2  +2CX6bxb2  cos  2y, 

+2 CX7bfb2  sin  2 /,  +2C]gb2b{  +2 Cx9b2bx  cos 2 yt  -2 C20b2b,  sin  2y,  -C2lY0u  sin  ye 

~YoubA  [  W  sin (ye  -  /, )  +  q2cb2  sin  (ye  +  y, )]}  (43c) 

y\  =  a,  +  {-^Ibl  sin /,  -bx  [C22ax  cosy,  +  2sxla>x  sin/,]  -2sna)fbx  sin/, 

‘\CDx(Q2OxQ1 

+b2  (CnA  cos/,  - 2 q20a?2  sin  /, )  +C246262  cos/,  +C256262  sin  /,  +  2C16b2b2 
+2b\  b2  (C27  cos  2/,  +  C28  sin  2/, )  +2C296236,  +2£&  (C30  cos  2/,  -  C31  sin  2/, ) 

+C32T0„  sin  ye  +  Y0ubxb2  [sxax  sin  (/,  -/,)-<72ft>2  sin  (ye  +  /,)]}  (43d) 

where  the  coefficients  C,  are  defined  in  Appendix  C.  Note  that  equations  (43)  contain  four 
unknowns  namely  the  response  amplitudes,  bx,  and  b2,  and  the  two  phase  angles,  ye  and  /, . 

Furthermore,  the  dependency  on  the  external  and  internal  detuning  parameters  explicitly  appears 
in  equations  (43c)  and  (43d).  The  solution  of  these  equations  is  obtained  numerically  for 
different  values  of  excitation  amplitude,  at  the  critical  and  post-critical  flutter  speeds.  In  both 
cases,  the  response  is  estimated  in  the  neighborhood  of  combination  parametric  resonance,  i.e., 
Q  » co,  +co2 ,  internal  frequency  ratio,  (0,  /co2  «  1 ,  external  excitation  detuning  parameter,  ct,  =  0, 

internal  resonance  detuning  parameter  cr ,  =0.061,  and  system  parameters  p  =  l/38,  a  =  -0.36, 
xa  =0.024,  ra  =  0.62 ,  r  =  0.5  =  0.001 ,  Qu  =  0.001 ,  </  =  0.1. 
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4.3.1  Response  at  Critical  Flutter  Speed 

In  the  absence  of  parametric  excitation,  the  solution  of  these  equations  reveals  the  coexistence  of 
different  fixed  points  implying  limit  cycle  oscillations  (LCO).  Each  fixed  point  is  created  from  a 
certain  domain  of  initial  conditions.  There  are  three  possible  values  of  the  fixed  point  and  each  is 
achieved  for  certain  values  of  initial  conditions.  The  three  fixed  points  are 

{w,a}  =  {0,0},  {0.15, -0.11},  and  {1.4,0.325} 

The  domains  of  attraction  of  each  of  the  above  fixed  points  are  shown  in  Figure  5(a),  where  the 
empty  space  is  belonging  to  the  zero  fixed  point. 

Under  parametric  excitation,  the  same  scenario  is  preserved  up  to  an  excitation  level 
Y0  =  0.00875  at  which  the  response  experiences  Hopf  bifurcation  around  each  fixed  point  in 

addition  to  a  new  fixed  point  that  also  experiences  Hopf  bifurcation.  At  T0  =  0.00875 ,  there  are 

four  different  values  of  initial  conditions,  which  lead  to  four  different  attractors.  These  initial 
conditions  are: 


{«0,a0,w0,a0}  =  {0.5,  0.04,  0, 0} ,  {0.5, 0.05, 0, 0},  {0.5, 0.7,  0,  0} ,  and  {0.45, 0.07,  0,  0} . 


The  domains  of  attraction  corresponding  to  this  case  are  shown  in  Figure  5(b). 

Over  another  range  of  excitation  amplitude,  0.00875  <  T0  <  0.0125  ,  the  response  experiences 
Hopf  bifurcation  for  relatively  high  range  of  initial  conditions  below  which  the  response  settles 
at  the  static  equilibrium  position.  For  example,  under  excitation  amplitude,  Y0  =0.01 ,  the  basin 

of  attraction  of  this  response  is  shown  in  Figure  5(c)  and  it  reveals  that  the  basin  of  initial 
conditions  that  leads  to  non-zero  response  is  shifted  away  to  either  higher  positive  values  or 
negative  values.  At  excitation  amplitude  T0=  0.0125  an  additional  bifurcation  around  a  new 
point  takes  place  depending  on  the  initial  conditions.  For  example,  the  two  sets  of  initial 
conditions  {0.3, -0.1,  0,  0}  and  {0.5, 0.1, 0, 0} ,  yield  two  different  periodic  attractors  under 

T0  =  0.014 .  Figure  5(d)  shows  the  basins  of  attraction.  In  the  first  and  third  quadrants  the  basins 

of  attraction  are  diminished  to  one  positive  or  negative  point,  i.e.,  {0.5, 0.1, 0, 0}  or 

{-0.5, -0.1,0, 0}.  These  two  sets  of  initial  conditions  lead  to  large  amplitude  oscillations. 

Another  two  basins  of  attraction  located  in  the  second  and  fourth  quadrants  lead  to  another 
attractor  of  very  small  amplitude. 

As  the  excitation  amplitude  increases,  additional  bifurcations  take  place  in  the  form  of  cascade  of 
period  doubling  depending  on  initial  conditions.  Figure  6(a)  shows  samples  of  time  history 
records,  phase  diagrams,  and  FFT  plots  under  excitation  amplitude  T0  =0.015,  and  initial 

conditions  {0.5,  -0.1, 0,  0} .  It  is  seen  that  the  bifurcation  takes  place  in  the  form  of  period 

doubling.  The  domain  of  attraction  of  this  attractor  is  slightly  enlarged  than  the  corresponding 
case  of  Figure  5(d).  There  is  another  periodic  attractor  of  period  two,  shown  in  Figure  6(b)  that 
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oscillates  with  larger  amplitude  than  the  first  attractor  and  emanates  from  one  point  of  initial 
conditions  in  the  first  and  third  quadrants  shown  in  Figure  5(e). 

Under  excitation  amplitude  Y0  =  0.016,  the  response  experiences  another  bifurcation  of  period 

four.  There  are  two  different  attractors  similar  to  those  displayed  in  Figure  5(e).  Any  further 
increase  of  excitation  amplitude  results  in  chaotic  motion  as  shown  in  Figure  7,  which  is 
generated  under  excitation  amplitude  T0  =0.017  and  different  domains  of  attraction.  The 

spectrum  of  the  time  history  records  is  continuous.  The  Poincare  maps  displayed  in  Figures  7(a) 
and  (b)  do  not  cover  a  smooth  closed  curve  but  constitute  a  set  of  randomly  scattered  points 
implying  that  the  motion  turned  out  to  be  chaotic.  The  Poincare  maps  are  obtained  based  on  the 
period  of  the  first  return.  The  domains  of  attraction  corresponding  to  the  small  amplitude 
oscillation  attractor  are  expanded  as  shown  in  Figure  5(f). 

Figure  8  presents  the  bifurcation  diagram  by  taking  the  excitation  amplitude  as  the  control 
parameter.  It  is  seen  that  for  a  region  of  initial  conditions,  the  parametric  excitation  acts  as  a 
stabilizer  source  of  the  wing  flutter  where  under  these  regions  the  response  achieves  its  static 
equilibrium  position.  Outside  this  domain  of  initial  conditions,  the  equilibrium  position  loses  its 
stability  and  the  response  may  possess  fixed  points,  Hopf  bifurcation,  cascade  of  period 
doubling,  and  eventually  chaotic  motion.  It  is  seen  that  at  each  bifurcation  point,  branches  of 
symmetric  and  asymmetric  periodic  solutions  meet  in  the  form  of  supercritical  symmetry¬ 
breaking  bifurcation. 

It  is  obvious  that  the  response  exhibits  multiple  attractors,  each  with  its  own  domain  of 
attraction.  One  or  two  of  these  attractors,  namely  the  zero  response  or  small  amplitude  LCO 
about  a  small  mean,  are  superior  and  desirable  in  the  post  flutter  region.  In  order  to  achieve  that 
desirable  performance  one  may  use  one  of  the  current  techniques  used  to  control  chaos  (see,  e.g., 
[53]).  For  the  present  case,  the  domains  of  attractions  of  undesirable  performance  can  be  shifted 
away  only  by  increasing  the  parametric  excitation  amplitude.  Figure  9  shows  the  dependence  of 
the  percentage  of  domains  of  attraction  that  lead  to  non-zero  flutter  oscillation  on  the  excitation 
amplitude.  It  is  seen  that  as  the  excitation  amplitude  increases  both  the  domains  of  attraction  for 
low  and  large  response  amplitude  decrease  until  critical  excitation  amplitude,  Ycr  »  0.0 1375, 
above  which  the  parametric  excitation  suppresses  the  wing  flutter. 

4.3.2  Response  at  Post-Critical  Flutter  Speed 

At  a  flow  speed  that  is  slightly  higher  than  the  critical  flutter  speed,  Ux/b(oa  =5.02,  the  wing 

enters  in  the  post  flutter  region.  In  the  absence  of  parametric  excitation  the  wing  experiences 
Hopf  bifurcation  of  different  amplitude  oscillations  depending  on  initial  conditions.  For  all 
values  of  initial  conditions  occupied  in  Figure  10(a),  there  are  two  domains  of  response 
amplitudes.  Low  oscillation  amplitudes  shown  by  the  region  of  small  solid  triangles,  ▲,  and 
large  amplitudes  shown  by  solid  squares,  ■.  Under  very  small  excitation  amplitude,  the  response 
experiences  multi-periods.  Over  the  excitation  amplitude  range,  0.001  <  Y0  <  0.003 ,  the  response 
experiences  multi-periods  with  growing  amplitudes.  For  example,  under  excitation  amplitude, 
T0  =0.0025  ,  Figures  1 1(a)  and  (b)  show  time  history  records,  phase  portraits,  Poincare  mapping 
and  FFT  plots  for  two  different  sets  of  boundary  conditions. 
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Over  another  excitation  amplitude  range,  0.003  <Y0<  0.005  ,  there  are  two  domains  of  attraction 

one  of  them  leads  to  small  amplitude  oscillations  with  multi-period  as  in  the  previous  range, 
while  the  other  leads  to  periodic  then  cascade  of  period  doubling  with  high  amplitude 
oscillations.  Figures  12(a)  and  (b)  show  samples  of  time  history  records,  phase  diagram,  and  FFT 
for  the  two  cases  under  excitation  amplitude  Y0=  0.004.  Figure  12(a)  includes  also  Poincare 
maps. 

The  next  region  of  excitation  amplitude,  0.005  <  Y0  <  0.01 ,  is  characterized  by  chaotic  motion  for 

all  possible  initial  conditions.  However,  one  set  of  initial  conditions  leads  to  a  chaotic  attractor 
with  small  amplitude  oscillations,  while  the  other  set  leads  to  large  amplitude  oscillations  as 
shown  in  Figures  13(a)  and  (b),  respectively  for  T0  =0.0075 .  There  is  a  small  window  around 

Y0  =  0.01  characterized  by  period  doubling  for  all  possible  initial  conditions. 

A  new  regime  of  response  behavior  emerges  at  excitation  level  range,  0.01  <  T0  <  0.014 

characterized  by  a  train  of  spikes  known  as  the  “firing'  state  [54]  for  certain  regions  of  initial 
conditions.  The  remaining  initial  conditions  lead  to  periodic  attractor.  For  example,  under 
excitation  amplitude,  Y0  =  0.013,  Figures  14(a)  and  (b)  show  typical  time  history  records  of  the 

firing  state.  Figures  14(c)  and  (d)  show  the  periodic  regime.  Figure  10(b)  show  the  domains  of 
attraction  of  these  four  attractors.  Note  the  period  between  the  spikes  in  the  firing  state  varies 
substantially.  As  the  excitation  amplitude  increases  the  period  between  the  spikes  (called  the 
“ refractory ”  or  “ recovery ”  period)  increases  and  for  the  some  initial  condition  the  recovery 
period  becomes  infinitely  large  declaring  a  stabilization  effect.  Figures  15(a)-(c)  show  samples 
of  time  history  records  belonging  to  these  regimes.  Figure  10(c)  shows  the  domains  of  attraction 
that  lead  to  these  different  attractors.  At  excitation  levels,  T0  >  0.015,  the  response  achieves  the 
zero  equilibrium  position. 

Figures  16(a)  and  (b)  show  the  bifurcation  diagram  which  summarizes  all  the  above  stated 
regimes.  In  order  to  appreciate  the  stabilization  effect  of  parametric  excitation,  Figure  17  gives 
the  dependence  of  the  percentage  of  the  area  of  domain  of  attraction  on  the  excitation  amplitude. 

5.  NUMERICAL  SIMULATIONS 

The  purpose  of  the  numerical  simulation  of  the  original  equations  of  motion  is  to  validate  the 
multiple  scales  predictions.  Equations  of  motion  (15)  and  (16)  are  numerically  integrated  for 
different  values  of  excitation  amplitude.  The  numerical  integration  is  carried  out  in  the  absence 
of  air  flow  and  in  the  presence  of  air  flow  at  the  critical  flutter  speed  for  p  =  l/38,  a  =  -0.36, 

*a  =0.024,  ra=  0.62,  r  =  0.5  ,C,a  =  0.001 ,  =0.001,  d  =  0.1.  The  following  subsections 

present  just  representative  samples  of  the  numerical  simulation  and  are  not  considered 
exhaustive. 

5.1  ZERO  FLOW  SPEED 

Under  parametric  excitation  with  frequency  close  to  the  sum  of  the  first  two  mode  frequencies, 
Q  «  a),  +  co2 ,  and  normal  modal  frequencies  ratio  to,  /co2  *  0.5 ,  the  wing  response  is  obtained  for 

different  excitation  amplitudes.  Under  all  values  of  excitation  amplitude,  the  response  of  both 
bending  and  torsion  exhibits  periodic  modulated  signals  as  shown  in  Figures  18(a)-(d).  At 
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relatively  low  excitation  amplitude,  the  modulation  effect  is  very  weak  as  demonstrated  in  Figure 
18(a)  and  becomes  significant  as  the  excitation  amplitude  increases.  The  simulations  are  carried 

out  for  initial  conditions  «(0)  =  0.5,  a(0)  =  0.1,  w'(0)  =  0,  a'(0)  =  0.The  bifurcation  diagram 
is  shown  in  Figure  19  which  is  obtained  for  initial  conditions  i7(0)  =  0.5,  a(0)  =  0.1,  w'(0)  =  0, 
a'(0)  =  0. 

5.2  CRITICAL  FLUTTER  SPEED 

Under  flutter  flow  speed  Ux  /b(oa  =  4.965  and  parametric  excitation  with  frequency  close  to  the 
sum  of  the  first  two  mode  frequencies,  Q  =  (ox  +  0)2 ,  and  internal  frequencies  having  the 
relationship  co2  =  to, +0.061 ,  the  wing  response  is  obtained  for  different  values  of  excitation 
amplitudes  and  initial  conditions.  Under  zero  excitation  amplitude,  and  depending  on  initial 
conditions,  there  exist  two  different  LCOs.  At  excitation  amplitude  T0  =  0.00875 ,  and  depending 

on  initial  conditions  the  response  either  preserves  the  zero  equilibrium,  or  possesses  multi-period 
oscillations  as  shown  in  Figure  20(a).  These  two  attractors  are  maintained  up  to  excitation 
amplitude  0.013  above  which  the  response  is  attracted  to  either  the  zero  response  or  to  a  chaotic 
motion  as  shown  in  Figure  20(b)  for  T0  =  0.014 .  Note  that  the  shown  Poincare  maps  are  obtained 

based  on  the  excitation  period.  As  the  excitation  amplitude  increases,  say  to  T0=  0.016,  the 
chaotic  regime  is  characterized  by  occasional  spikes  followed  by  relaxation  period  as  shown  in 
Figure  20(c).  The  zero  response  is  not  shown  for  this  case.  At  excitation  amplitude  T0  =0.017 

the  non-zero  response  becomes  more  ordered  after  experiencing  spiky  transition  period  as  shown 
in  Figure  20(d). 

5.3  POST-CRITICAL  SPEED 

At  slightly  higher  value  of  flow  speed  than  the  critical  flutter  speed,  say,  U00/b(oa  =5.02,  the 

numerical  simulation  is  presented  for  selected  values  of  excitation  amplitude.  At  zero  excitation 
amplitude  the  wing  wings  experiences  LCOs  in  bending  and  torsion.  The  zero  attractor 
disappears  for  all  excitation  levels  up  to  excitation  amplitude  of  0.013.  As  the  excitation  level 
increases  from  a  low  value  up  to  0.013,  the  response  time  history  records  experience  multi- 
periods  and  Figures  21(a)-(d)  show  the  sequence  of  development  of  the  response  time  history 
records  for  different  values  of  excitation  amplitude.  For  excitation  levels  of  0.014,  and  above,  the 
zero  attractor  emerges  for  certain  domain  of  attraction  in  addition  to  a  multi-period  attractor 
depending  on  the  initial  conditions.  For  example.  Figures  22(a)  and  (b)  show  time  history  records 
at  excitation  level  0.01 5  for  two  different  sets  of  initial  conditions.  The  domains  of  attraction  that 
lead  to  zero  attractor  are  enlarged  as  the  excitation  amplitude  increases  and  this  confirms  the 
results  predicted  by  the  multiple  scales  method. 

5.  CONCLUSIONS 

The  nonlinear  flutter  of  a  cantilever  wing  is  studied  analytically  and  numerically  in  the  absence 
and  presence  of  parametric  excitation.  In  the  absence  of  nonlinearities,  the  regions  of  parametric 
instability  are  obtained  for  different  values  of  flow  speed. 

•  At  the  critical  flutter  speed  the  bottom  of  instability  region  touches  the  frequency  axis. 
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•  Below  and  above  the  critical  speed  the  instability  regions  move  away  from  the  frequency 
axis. 

•  As  the  excitation  amplitude  increases  the  responses  experiences  a  cascade  of  period 
doubling  and  eventually  chaotic  motion. 

•  At  critical  excitation  amplitude  the  response  exhibits  stabilization  state  to  the  zero 
equilibrium  attractor  over  a  wide  range  of  initial  conditions.  This  new  and  unexpected 
feature  motivated  the  authors  to  add  a  review  assessment  of  parametric  excitation- 
induced  stability. 

•  The  stabilization  effect  was  also  manifested  at  a  flow  speed  that  is  slightly  higher  than  the 
critical  flutter  speed.  However,  the  stabilization  effect  was  preceded  by  new  phenomena 
such  as  firing  and  recovery  states. 

•  The  numerical  simulation  of  the  original  equations  of  motion  has  confirmed  the  multiple 
scales  findings. 


Appendix  A:  Elements  of  equation  (22) 
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AA(ka ,  to)  =  co2  (  Ar  +  iA, )  +  2/r^tn  +  r2, 

DD(ka,G))  =  o)2  (DR  +  iDj) , 

where 
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Appendix  B:  The  parameters  p(,  i  =  1,2, 3, 4  in  equations  (27)  and  (28) 
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Appendix  C:  Coefficients  C,  in  equations  (43) 

^1  =  3<7|5  +  (7ll®2  >  ^2  =  ^7®!  ■*"9l8®2  » 


C3=9,|6+3^Wl  -?6<*V°I 
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4  *7|9  +  —  ’  Cj  —  _  +  <7|q(02  _  ,  C6  —  3S|4  +  •S’4C0,  ,  C7  —  S|9®|  +  SgC02 

CO,  (0,  CO, 

Cs  =  s3+  3s5co2  -  s10w,co2  ,  C9  =  -s,  +  sit  ,  C10  =  +  56co,  -  5,  ^ 

(0,  (02  C02 

-  6-5,4  +2-s4<y,  ,  C,2  =69,5  +  2quco2 , 

C\  3=2  (353fij,  +  39l6«2  +  (?9  +  s10 )  +  (q6+ss)  (dxco\  ) 

C|4  —  2[3.y,9<y,  +  (<y7  — ^8)<y|<y2  — 3^|g<y2  J,  C|5  =  ^2s^0)^  ■+■  2^9co,  +  3qt 4<y2  +  <74*^1  g>2  j 

^16  =  ^16^1  Sgty  +  s6(Ot  o2 ,  Cl7  =  57<y,  —  ,  Clg  =  2<73<y2  +  2<75<y2  +  3s1J<y|  +  Sn<y2<y, 

^19  —  #3^2  —  ^5^2  "*"  *7l0^2  *  C2fl  =  ^g<y2  _  ^19^1^2  >  ^21  =  S2®\^\  +  *71^2^2  »  ^22  =  ^14  "*"  2^4^|2  , 

C23  =  6?15  +  2tfl 1^22  ’  C24  =  2  [-3-Sjfi>,  +  3^,6S2  +  (g,  -  5|0  )  6>,26>2  +  (ft  -  i5  )  fi>,6>2  ] 

C25  =  2  [-3S|96>|2  +  (<jr7  +  5g )  ft>,ft>2  - 3qna>]  ] ,  C26  =  [-2s16g>,  -  2s9tf5,3  +  3qHco2  +  q^a)2  ] 

C21  =  “SI6^1  +  'W  -  S6®1®2  »  C28  =  -'W  +  *1*^2  *  C29  =  2 ^3 ^2  +  2ftfi>2  ~  3s,5tf>,  ~  •*,  * 

C30=^3^2_^5^2  +9l0^22^n  Ql  =  “  9|9^2  »  Cjj  =  ~  q&b]  . 
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Figure  1 .  Schematic  diagram  of  the  analytical  model,  coordinates  axes,  and  deformations. 

(a)  Plate-like  cantilever  wing  showing  the  coordinate  frame 

(b)  Locations  of  different  axes  and  deformations  in  bending,  torsion,  and 
longitudinal. 
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Figure  2.  Dependence  of  eigenvalue  components  on  the  flow  speed  parameter  showing:  (a)  the 
influence  of  the  imaginary  component  of  circulation  function,  G  *  0  and  (b)  its 
absence,  G  =  0 . 

Eigenvalue  components  in  the  presence  of  aerodynamic  and  structure  damping 
.  Eigenvalue  components  in  the  absence  of  aerodynamic  and  structure  damping. 


Figure  3.  Stability  boundaries  in  the  neighborhood  of  parametric  combination  resonance  at  zero 
flow  speed  and  at  critical  flow  speed. 
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Figure  4.  Bifurcation  diagram  at  zero  critical  flow  speed  estimated  at  Q  =  oj,+co2,  co,  =  0.5w2, 
for  cre  =0,  i7(0)  =  0.1 ,  a(0)  =  0.1,  w'(0)  =  0,  a’(0)  =  0. 

(a)  Dependence  of  bending  response  on  excitation  amplitude. 

(b)  Dependence  of  torsion  response  on  excitation  amplitude. 
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Figure  5.  Domains  of  attractions  for  different  values  of  excitation  amplitude. 

(a)  ro=0;  (b)  Y0  =  0.00875;  (c)  K0  =  0.01 ;  (d)  Y0=  0.014;  (e)  Y0  =  0.015;  (0 
Y0  =0.017 
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Figure  6(a) 
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Figure  6(b). 

Figure  6.  Time  history  records,  phase  portraits,  FFT  plots  at  critical  flow  speed 
£4, / =  4.965 ,  for  Y0  =  0.015 ;  Q  =  co,  +to2 ,  co,  =  to2 ,  =  0,  cr,  =  0.061 ,  and  for 

two  sets  of  initial  conditions. 

(a)  «(0)  =  0.5,  <x(0)  =  -0.1,  w’(0)  =  0,a'(0)  =  0; 

(b)  w(0)  =  0.5,  a(0)  =  0.1 ,  u  '(0)  =  0,a'(0)  =  0 


126 


-0-1  o  0.1  0.2  0.3  0.4  0.5 


a 


a 


-0.2  -0.15  -0.1  -0.05 


0.2  0.4  0.6  0.8  1 


co/co„ 


C0/(0« 


tf 


Figure  7(a) 
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Figure  7(b). 

Figure  7.  Time  history  records,  phase  portraits,  Poincare  maps,  and  FFT  plots  at  critical  flow 
speed  Um  /  b(Oa  =4.965  ,  excitation  amplitude  Y0  =  0.017  ;  Q  =  {0,+{02,  to,=co2, 

oe  =  0,  ct,  =  0.061 ,  for  two  different  sets  of  initial  conditions 

(a)  m  (0)  =  0.5 ,  a(0)  =  -0.1,  u' '(0)  =  0,a’(0)  =  0 

(b)  ?7(0)  =  0.5,  a(0)  =  0.1,  M'(0)  =  0,a'(0)  =  0 
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Figure  8.  Bifurcation  diagrams  at  critical  flow  speed  Ux/bc oa  =4.965,  Q  =  co,+o)2,  to,  =co2, 
<7^=0,  ct(  =  0.061 . 

(a)  Dependence  of  bending  amplitude  on  excitation  amplitude. 

(b)  Dependence  of  torsion  amplitude  on  excitation  amplitude. 
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Figure  9.  Stabilization  effect  of  parametric  excitation  at  critical  flow  speed  Ux/b(oa  =4.965 

showing  the  percentage  of  domains  of  attraction  for  large  amplitude  response  (solid 
diamond  curve)  and  large  plus  small  amplitude  responses  (dotted  solid  circle  curve). 
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Figure  10.  Domains  of  attraction  at  post-critical  flow  speed  Ux/b(oa  =5.02  for  f2  =  (o,+co2, 
co1=co2,  CTe=0,  ct,  =0.061,  and  different  values  of  excitation  amplitude,  (a) 
ro=0;(b)  Y0  =0.013;  (c)  Y0  =0.014. 
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Figure  1 1(b). 

Figure  11.  Time  history  records,  phase  portraits,  Poincare  maps  and  FFT  plots  at  post-critical 
speed  Ua/b(oa=  5.02,  Y0  =0.0025;  ft  =  co,+to2,  co,=co2,  =0,  a,  =0.061,  and 


for  two  different  sets  of  initial  conditions. 

(a)  i7 (0)  =  0.2,  a(0)  =  0.01 ,  M’(0)  =  0,a'(0)  =  0 

(b)  m(0)  =  0.5,  a(0)  =  0.1,  «'(0)  =  0,a'(0)  =  0 
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Figure  12(a) 
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Figure  12(b). 

Figure  12.  Time  history  records,  phase  portraits,  Poincare  maps  and  FFT  plots  at  post-critical 
speed  Un/ba oa=5.02,  T0  =0.004;  Q  =  co,+co2,  co,=(o2,  crc=0,  a, =0.061,  and 
for  two  different  sets  of  initial  conditions. 

(a)  m(0)  =  0.2,  a(0)  =  0.01 ,  w’(0)  =  0,a'(0)  =  0 

(b)  w(0)  =  0.5,  a(0)  =  0.1 ,  i7'(0)  =  0,a'(0)  =  0 
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Figure  13(a). 


135 


a 


a 


0.05 

0.025 

0 

-0.025 

-0.05 

-0.075 

-0.1 


Figure  13(b). 

Figure  13.  Time  history  records,  phase  portraits,  Poincare  maps  and  FFT  plots  at  post-critical 
speed  U^/b^  =5.02,  Y0  =0.0075;  Q  =  co1+co2,  co,=co2,  a,  =0,  a, =0.061,  and 


for  two  different  sets  of  initial  conditions. 

(a)  w  (0)  =  0.2,  a(0)  =  0.01 ,  w,(0)  =  0,a,(0)  =  0 

(b)  w (0)  =  0.5 ,  a(0)  =  0.1 ,  w '(0)  =  0,a’(0)  =  0 
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Figure  14.  Time  history  records  showing  the  influence  of  initial  conditions  on  generating  firing 
(a,  b)  and  periodic  (c,  d)  states  at  post-critical  speed  UK/b(oa  =5.02,  T0  =0.013; 

Q  =  to,  +  co2 ,  to,  =  co2 ,  a,  =  0 ,  a,  =  0.06 1 . 

(a)  i7(0)  =  0.1,  a(0)  =  0.02,  «’(0)  =  0,a'(0)  =  0 

(b)  m(0)  =  0.  1,  a(0)  =  0.02,  «'(0)  =  0,a'(0)  =  0 

(c)  m(0)  =  0.2,  a(0)  =  0.02,  t7'(0)  =  0,a’(0)  =  0 

(d)  m(0)  =  0.4,  a(0)  =  0.1 ,  i7'(0)  =  0,a'(0)  =  0 
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Figure  15.  Time  history  records  showing  the  influence  of  initial  conditions  on  generating  firing 
(a,  b)  and  periodic  (c,  d)  states  at  post-critical  speed  Ua0/bcoa  =5.02,  Y0  =0.014; 

Q  =  co,  +cd2  ,  co,  =  co2 ,  ne  =  0,  a  =  0.061 . 

(a)  *7(0)  =  0.3,  a(0)  =  0.02 ,  *7’(0)  =  0,a’(0)  =  0 

(b)  w(0)  =  0.35 ,  a(0)  =  0.02 ,  w'(0)  =  0,a'(0)  =  0 

(c)  *7(0)  =  0.25,  a(0)  =  0.07 ,  i7'(0)  =  0,a’(0)  =  0 
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Figure  16.  Bifurcation  diagrams  showing  different  regions  of  response  behavior  at  post-critical 
flow  speed  UK /b(oa  =  5.02 ,  Q  =  co,  +co2 ,  co,  =  co2 ,  ae  =  0,  a, .=  0.061 

(a)  Dependence  of  bending  amplitude  on  excitation  amplitude 

(b)  Dependence  of  torsion  amplitude  on  excitation  amplitude 
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Figure  17.  Stabilization  effect  of  parametric  excitation  at  post-critical  flow  speed 
Um  /  be oa  =  5.02  showing  the  percentage  of  domains  of  attraction  for  large  amplitude 

response  (diamond-solid  curve)  and  large  plus  small  amplitude  responses  (star- 
dotted  curve). 
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Figure  18.  Response  time  history  records  according  to  numerical  simulation  at  zero  critical  flow 
speed  for  different  values  of  excitation  amplitude  and  for  the  same  initial  conditions 


m (0)  =  0.5 ,  a(0)  =  0.1 ,  i7'(0)  =  0,  a’(0)  =  0. 


(a)  Y0  =  0.001 ;  (b)  Y0  =  0.005 ;  (c)  Y0  =  0.01 ;  (d)  Y0=  0.01 5 . 
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Figure  19.  Bifurcation  diagram  according  to  the  numerical  simulation  of  equations  of  motion  at 
zero  flow  speed  and  fi  =  co,+(o2,  (o,=0.5co2,  «(0)  =  0.5,  a(0)  =  0.1,  i7’(0)  =  0, 

a'(0)  =  0 . 

(a)  Dependence  of  bending  response  amplitude  on  excitation  amplitude 

(b)  Dependence  of  torsion  response  amplitude  on  excitation  amplitude 
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Figure  20(a) 
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Figure  20(b) 
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Figure  20(c) 
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Figure  20(d). 

Figure  20.  Selected  samples  of  numerical  simulation  time  history  records,  phase  portraits, 
Poincare  maps  (for  b,c)  and  FFT  plots  for  different  values  of  excitation  amplitude  and 
different  sets  of  initial  conditions,  at  critical  flow  speed. 

(a)  Y0  =0.00875,  w(0)  =  0.5,  cc(0)  =  0.04,  *7'(0)  =  0,  a’(0)  =  0 

(b)  Y0  =0.014,  *7(0)  =  0.5,  a(0)  =  0.1,  w'(0)  =  0,  a'(0)  =  0 

(c)  Y0  =  0.016,  m(0)  =  0.6,  a(0)  =  0.1,  i7'(0)  =  0,  a'(0)  =  0 

(d)  Y0  =0.017,  i7(0)  =  0.6,  a(0)  =  0.1 ,  «'(0)  =  0,  a’(0)  =  0 
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Figure  21.  Sample  of  numerical  simulation  time  history  records  at  pos-critical  speed 
Un/b(da=  5.02,  fi  =  co,  +  co2 ,  co2  =  co,  +0.061 ,  i7(0)  =  0.1,  ct(0)  =  0.01,  w'(0)  =  0, 

a'(0)  =  0 ,  and  for  different  values  of  excitation  amplitude. 

(a)  Y0  =  0.001 ;  (b)  Y0  =  0.005 ;  (c)  Y0  =  0.0075 ;  (d)  Y0  =  0.01 


147 


Figure  22(a). 
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Figure  22(b). 

Figure  22.  Samples  of  numerical  simulation  time  history  records  at  post-critical  flow  speed 
Ua /b(oa  =  5.02 ,  excitation  amplitude  Y0  =0.015  for 

(a)  a  set  of  initial  conditions  that  lead  to  zero  response:  i7(0)  =  0.2,  a(0)  =  0.01, 
w’(0)  =  0,  a'(0)  =  0 

(b)  a  set  that  lead  to  quasi-periodic  response:  «(0)  =  0.3,  a(0)  =  0.01,  w"(0)  =  0, 
a'(0)  =  0 
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